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Abstract. Recently, we studied the technologically important problem of pe¬ 
riodically forced spheroids in.simple shear flow and demonstrated the existence 
of chaotic parametric regimes. Our results indicated a strong dependence of the 
solutions obtained on the aspect ratio of the spheroids, which can be used to sep¬ 
arate particles from a suspension. In this paper we demonstrate that controlling 
the chaotic dynamics of periodically forced particles by a suitably engineered 
novel control technique, which needs little information about the system and 
is easy to implement, leads to the possibility of better separation. Utilizing the 
flexibility of controlling chaotic dynamics in a desired orbit irrespective of ini¬ 
tial state, we show that it is theoretically possible to separate particles much 
more efficiently than otherwise from a suspension of particles having different 
shapes but similar sizes especially for particles of aspect ratio > 1.0. The 
strong dependence of the controlled orbit on the aspect ratio of the particles 
may have many applications such as in the development of computer-controlled 
intelligent rheology. The results suggest that control of chaos as discussed in 
this work may also have many applications. 
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1. Introduction 

There has been appreciable interest in the dynamics of suspensions of small dipolar par¬ 
ticles stabilized by surface active agents in a variety of linear flows under the effect of 
alternating or rotating external fields. This problem has been studied by a number of au¬ 
thors and a number of practical applications such as magnetofluidization (Buevich et al 
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1984), magnetostriction of ferromagnetic particle suspensions (Ignatenko et al 1984), the 
growth of single crystals from a melt by Czochralski’s method (Zibol'd et al 1986), rhe¬ 
ological properties of feiTomagnetic colloids (Tsebers 1986) and the characterization of 
magnetorheological suspensions (Cebers 1993; Petrikevich & Raikher 1984; Kashevskii 
1986; ShuPman et al 1986) have been discussed. Some of the above authors consider ro¬ 
tatory external fields which can be taken as a superposition of two orthogonal alternating 
fields. As can be seen from the relevant literature cited above, the problem discussed in 
this paper is of interest in a number of technological applications. For a brief literature 
review of the importance, relevance and novelty of this problem, the reader is referred to 
our earlier paper (Kumar et al 1995) and the references therein. In a preliminary analysis of 
this problem we have demonstrated (Kumar et al 1995) the existence of chaotic dynamics 
and a strong dependence on the aspect ratio of the particles on the solution of the problem 
in the chaotic regime. This strong dependence was considered to demonstrate the potential 
of this problem to develop techniques for effective separation of particles by aspect ratio 
to yield well-characterized suspensions for testing theories of the above processes and for 
developing suitable theories for suspensions of particles of various aspect ratios and sizes. 
It is also of interest to study the possibility of controlling the chaos present in the dy¬ 
namics of periodically forced particles with a view to develop potential applications such 
as control of magnetorheological suspensions, magnetostriction of ferromagnetic particle 
suspensions, magnetodielectric material suspensions, ferromagnetic colloidal suspensions 
and the possibility of computer-controlled intelligent rheology. 

In order to separate out the effects of particle size and shape on the properties of the 
suspension, it would be desirable to obtain suspensions of particles having a narrow dis¬ 
tribution of shapes and sizes. A narrow distribution of sizes can be obtained through 
appropriate separation techniques such as filtration. However, we are not aware of any 
effective technique to separate particles having the same size but different shapes. 

There has been considerable interest recently in the use of chaos and controlling chaos 
as a tool to develop new processing technologies (Zumbmnnen et al 1996). Practically 
speaking, chaos may be desirable or undesirable depending on the application desired. 
In many practical situations, in order to improve system performance, the chaotic system 
must be controlled to a periodic orbit or to a steady state. Thus the technique of controlling 
chaotic systems has received increased attention in the recent literature. In a physical 
apparatus, one can imagine that the system dynamics is to be changed in some way so that 
improved performance is obtained. One can achieve this by making appropriate changes in 
file system to achieve the desired objective. Similar or otherwise unattainable objectives can 
be attained by operating the system in a chaotic regime and then applying suitable control 
of chaos techniques. Since a chaotic attractor is a closure of unstable periodic orbits of 
different periods, any of these orbits can be stabilized to attain the desired objective by 
suitable control algorithms. 

The existence of chaotic parametric regimes in the problem considered in this work and 
the possibility of control of its chaotic behaviour (Ott et al 1990) allows the possibility of 
many applications including the development of computer controlled intelligent rheology. 
Another possible application of controlling the chaos present in this problem is further 
fine-tuning of the separation of the particles by aspect ratio than is possible with other 
methods. Since the chaotic phase space is an unlimited reservoir of periodic solutions 


(Ott et al 1990), a particle of a particular aspect ratio can be forced to oscillate in any 
one of the orbits and can be picked up from there effectively. Though the results presented 
herein represent a preliminary analysis of the problem considered in dilute suspensions, the 
existence of chaos and hence its controllability and the strong dependence of the results on 
the aspect ratio of the controlled orbit suggest that these results may have many potential 
applications. In this paper, as an example, we focus on the more efficient separation of 
particles by aspect ratio through control of chaos. 

2. Theory 

The general equation governing the dynamics of suspensions of small dipolar spheroids 
of finite aspect ratio in simple shear flow under the effect of externally induced periodic 
fields, as developed by Kumar et al (1995), leads to 
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U 2 = Pui(l — 2 u2^) — Qu\ + Ra 2 cos,(cot), 
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with ^ depending on the axis ratio being > 1 or < 1 , u is the rate of change of the 
orientation vector; is the particle axis ratio of the spheroid defined as = a/b, where 
a = polar radius and b = equatorial radius; u\,u2, M 3 denote the x, y, z components 
of the unit vector u along the axis of symmetry which determines the orientation of the 
spheroids, respectively; k\ , k2 and k^ are the x, y, z components of the external force k; 
0 ) is the frequency of the driver and t is the time. 

As can be seen, the equations depend on the aspect ratio of the particle and also depend 
on the magnitude and frequency of the external force field. The maximum scaled value of 
0 ), corresponding to the limit where the quasistatic Stokes’ equations are valid corresponds 
to J = 2Tc{re re~^). The three equations listed above can be reduced to a system of 
two nonlinear nonautonomous ordinary differential equations (ODEs), since u\, U2 and 
M 3 evolve with time, from the given initial condition |u| = 1 , the magnitude of the vector 
u remains at unity. Thus, (1) can be converted into equations for 6 and 0 in spherical 
co-ordinates as given by them, where 6 is the azimuthal angle and 0 is the polar angle 
corresponding to a given vector u as 

de P 

— = — sin 20 sin 20 
dt 2 

R[k] COS0 COS0 -k k 2 cos 0 sin0 — /C 3 sin 6] cos{cot). 
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We note that the equation for 0 has a singularity at 0 =0. The singularity of 0 can be 
removed by taking 0i = 0 sin 0. The equations for 0 and 0 when integrated numerically, 
did not lead to any problems with the singularity at 0 = 0 and we obtained the same results 
as that obtained by the other set of equations. For a more detailed discussion of the method 
of solution of (2) the reader is referred to Kumar et al (1995). In this paper we suggest, 
as an example of the importance of chaos control, a technique to separate particles having 
the same size but different shapes much more efficiently than was demonstrated by them. 

The fact that a chaotic solution eliminates the possibility of long-term prediction of 
system behaviour induced many reports in the literature of either quenching chaos or 
controlling chaos (Ott etal 1990; Ditto etal 1995). Since chaotic attractors have embedded 
within them a dense set of unstable periodic orbits, any one of the unstable periodic orbits 
can be stabilized to obtain otherwise unattainable system behaviour. The essential idea is 
that a chaotic system explores a relatively large region of state space and the system can be 
brought to a desired stable state to improve the performance of the separation technique by 
a suitable control algorithm. The first method (OGY - acronym formed from the authors’ 
names) of control of chaos proposed by Ott etal (1 990) generated appreciable interest in the 
literature of chaos. Thereafter, a large number of algorithms for controlling chaos have been 
reported in the literature (Christini & Collins 1995; Rhode et al 1995). Broadly speaking 
there are two classes of algorithms for controlling chaos, namely (1) Feedback methods 
(Ott et al 1990) and (2) Non-feedback methods (Giiemez et al 1994). The first method needs 
appreciable information about system behaviour but is comparatively simple to implement 
experimentally. The second method may be difficult to implement experimentally. 

The algorithm suggested in this work needs very little information about system be¬ 
haviour but is rather easy to implement experimentally. This algorithm can be classified as 
a non-feedback control method. Rajasekar & Lakshmanan (1992, 1993) have investigated 
the applicability and effectiveness of various approaches of controlling chaos in the BVP 
oscillator. Suppression of chaos by periodic parametric perturbations in an experimental 
set-up of a Duffing oscillator is also reported in the literature by Fronzoni et al {1991). la. 
the method proposed in this paper, the control parameter is perturbed for one period at fixed 
intervals of every integral multiple of the fundamental period. In the methods considered 
earlier in the literature (Lakshmanan & Murali 1996), the parameter was perturbed con¬ 
tinuously rather than at fixed intervals. To implement the control strategy reported in this 
paper, we apply a constant external force in addition to the periodic force for the duration 
of one period after every n — \ periods as shown in figure 1, where length of one period 
(T) is equal to length of the fundamental period of the periodic force. If the system is 
perturbed after n-l periods by applying an additional constant force, the system is found 
to be stabilized in a periodic orbit with a period equal to nT or with a period equal to an 
integral multiple oinT. 

For example, if a represents the magnitude of the additional constant force, then the 
control can be implemented by setting a = 0 if j #0 mod(n) and a 0 if j = 0 
mod(n) where j represents iteration number and the time interval between two successive 
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Figure 1. The schematic representation of the control strategy presented in this paper, 
where control is active during the period corresponding to the shaded area and not active 
during the period corresponding to the blank area. 


srations is equal to the fundamental period T. That is, the system evolves without any 
Lodification for periods which are not multiples of n and with modification for periods 
hich are multiples of «. A schematic representation of the control strategy is given in 
gure 1. The chioce of n depends on the period-m solution we want to stabilise. The integer 
can be either m or a divisor of m depending on the value of a and the choice of n. 

To incorporate the above idea of controlling the dynamics of the system, we modify 
1 ) governing the dynamics of the system by introducing an additional constant force in 
Idition to the periodic force along the direction of the periodic force. Let k \, ^3 

2 the X, y, z components of the additional constant force, k'. After scaling all quantities 
spearing in the equations similarly to Kumar et al (1995), (2) can be written as 

d0 P 

— = Y sin 26 sin If + i?[(cos 6 cosf k\ + cos 0 sin f k 2 — sin 0 kf) 

X cos(£yO] + R[(cos 6 cos f k\ + cos 9 sin 0 /c 2 — sin 6 k'^)], 
df R 

— — P cos 20 — <2 H -[(“ sin f k] + cos 0 k 2 )cos(cot) 

dr sin 9 

4 - (— sin 0 k[ + cos 0 k' 2 )] 

\0\, | 0 l< 7 r/ 2 . (3) 

It is noticed that the modified system of equations reduces to (2) when k'^ — k 2 = k'^ 
: 0 . 

The application of an additional external forcing (constant) to control a chaotic system 
solving under the effect of external forcing (periodic) has been reported in the literature 
..akshmanan & Murali 1996). In this method the constant force is applied continuously, 
►ur control strategy involves the application of the constant force for a period of finite 
mgth, T after a period of length, {n — l)T (where T = 27r/cu is the fundamental period 
f the periodic forcing). Then, we lift the control for a period of finite length, (n — l)T 
nd thereafter we apply the constant force of the same magnitude for a period of length, 

’ and again we lift the control for a period of finite length, {n — l)T. This process is 
jpeated. Under this process, the system is allowed to evolve according to the system of 
i) upto the (n — l)th forcing period and evolve according to the system of (3) at the nth 
arcing period. This process is repeated every nth period. While solving (3), this idea can 
e implemented by setting k[ = = ^3 = 0 if j # 0 mod(n) and k\ = ^ 0 

' 7=0 mod(n), where j represents iteration number and T =2n/co is the time interval 
etween two successive iterations. We note that the above equations for 9 and 0 decouple 
1 the absence of ki , k[ and k 2 ,k 2 - Hence the presence of an external force field with either 

1 or ic ripppccsirv tr* nhtjiin rhantir cnliitinns in this svstpm Tn oiir r.alr.nlatinns wft kf.nt 



stabilized to a penodic orbit or penoa equal to an integral multiple or n. un me omer nana, 
if the control is applied continuously, we lose the flexibility of conti'olling the system to 
an orbit of desired period. 


3. Results and discussion 

Recently Kumar et al (1995) have analysed the properties of the system without control. 
They tentatively identified chaotic regimes of the parameter k 2 keeping k\ = ko, = 0. 
As a first step in analysing the properties of the equations derived in this paper, we set 
ki = k^ = Q,k\ = k 2 — kj =0 and varied k 2 for particles of different aspect ratio within 
the range of ranging from 0.2 to 2.0 in steps of 0.2 and kept co equal to J. We expect 
the greatest complexity of the solutions of these equations for this particular choice of 
parameters, since ^2 is responsible for the greatest opposition to the hydrodynamic torque 
due to the imposed shear flow field. At ki = 0, Jeffery’s (1922) results are reproduced and 
all solutions of the equations starting from different initial conditions tend towards a fixed 
point in the stroboscopic plot. 

We selected ^2 = 12.0 and computed the solution of the equations for different 
ranging within 0.2 to 2.0 in steps of 0.2. For ^2 = 12.0 the system given by (2) behaved 
chaotically for all the considered except for equal to 1.8 and 2.0. The Lyapunov 
exponents of the attractors were evaluated and found to be positive. The time series and 
attractor of a typical trajectory are shown in figures 2a and b for the case of = 1.2. 

It was reported by Kumar et al (1995) that the results of their computations are very 
sensitive to the aspect ratio of the particle in some parameter regimes. In the case of 
constant external fields in the same parametric regimes they obtained regular behaviour. 
For re > 1.0, they obtained nearly the same fixed point for all initial conditions in the 
case of a constant force field. This indicates that in the sample application considered 
in this work, a periodic force field is necessary to effect particle separation for particles 
with re > 1.0. Tables 1, 2 and 3 give a sample of the results obtained for a periodic 
force of amplitude k 2 = 12.0, a constant force of amplitude k 2 = 12.0 and for zero 
force respectively. In this case independent separation of particles is possible only for 
particles of aspect ratio re equal to 1.2. The existence of chaotic dynamics in this system 
allows control of its dynamics to a desired orbit and thus suggests the possibility of better 
separation of particles for almost all particles of aspect ratio ranging within 0.2 to 2.0. 
This is difficult to obtain in the case of regular behaviour or chaotic behaviour without 
control. 

We stabilized the chaotic dynamics obtained at ^2 = 12.0 to a periodic behaviour of 
period equal to a multiple of the fundamental period, liz/coXo demonstrate the applicability 
of the method proposed. The time series of U 2 of the Poincare section (stroboscopic plot) 
computed from (2) for re = 1.6 is given in figure 3. It is demonstrated that the system can 
he controlled to periodic behaviour of any desired period by applying the same constant 
force. For example, the system could be controlled to period-2, period-3, period-4 and 
period-5 orbits by applying the same constant force equal to = 5 as can be seen from 
figures 4, 5, 6 and 7. This is obtained by setting k\ = = kj = 0 if 7 ^ 0 mod(/i) 
and k[ = k’ 2 ^ — 0, kl^ = 5.0 if 7 =0 mod(n) in the evolution of (3) when the system is 
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Table 1. Distribution of evolution of initially uniformly distributed particles of different aspect ratios 
for the case k 2 = 12; w = /; 5 < m < 49 and 100 < / < 1000, where m is the total number of particles 
in the grid on the average and I is the total number of occurrences of the grid. 
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to be stabilized in a period-n orbit. The system evolves according to the evolution (2) for 
y ^ 0 mod(n) and according to the evolution (3) for 7=0 mod(n). It is also possible 
to stabilize a given chaotic system to different n-period orbits by suitably changing the 
control parameter As an example, the time series of Poincare sections corresponding 
to two different 2-period solutions embedded in a chaotic system are given in figure 8. 

The control technique is applied after 2000 fundamental periods and the system is 
followed upto 7000 periods with control. Once control is applied the system stabilizes 
rapidly to the appropriate periodic orbit. The magnitude of the perturbation required to 
stabilize the system was small compared to the magnitude of the periodic force for all 
aspect ratios except — 0.2, 0.4 and 0.6. The magnitude of the constant force required 
for control depends on the aspect ratio of the particle and the desired period, n for aspect 
ratios less than or equal to 0.6. 

One important advantage of the control algorithm outlined in this work is the possibility 
of switching over to different periodic solutions during a given run. This implies that a 
system in chaotic dynamics can be stabilized to one particular periodic orbit for a given 
time and to another periodic orbit of entirely different period after a given time as shown in 
figure 9. In this example the control parameter = 5.0 is applied between 2001 and 3000 
periods to stabilize the system to a period-2 orbit. Then the control is lifted between 3001 
and 4000 periods and again the control parameter = 5.0 is applied between 4001 and 
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5000 periods to stabilize the system to a period-3 orbit. The control was removed again 
between 5001 and 6000 periods. Thus the system oscillates in a period-2 orbit between 
2001-3000 periods and then oscillates in a period-3 orbit between 4001-5000 periods as 
can be seen from figure 9. This figure also reveals the fact that once the control is lifted, the 
system returns to a chaotic system. Another advantage which is important from the point 
of view of the sample application proposed in this work is the possibility of changing the 
periodic behaviour to another orbit of the same period by suitably changing the magnitude 
of the applied constant force. This allows us to bring a particle having a definite aspect 
ratio to a desired orbit by changing control parameters. It is also noted that all initially 
uniformly distributed particles of a given aspect ratio within the range 0.8 < < 2.0 can 

be concentrated in a given grid by applying a periodic force with control as can be seen 
from tables 4, 5 and 6. 

The observations made in this work suggest the possibility of separating particles more 
efficiently based on control of chaotic dynamics. We applied a control force (constant 
force) of magnitude 1 < k 2 < 5. For this range of k!^, the chaotic dynamics of particles 
of all aspect ratio except = 0,2,0.4 and 0.6 could be controlled to a desired orbit. This 
range of control force seems to be sufficient for efficient separation of particles by shape 
since particles of these aspect ratios alone can be brought to a desired position. To develop 
quantitative results based on this observation, we divided the range of possible orientations 





















Table 3. Distribution of evolution of initially uniformly distributed particles of different aspect ratios 
for the case ^:2 = 0; o) = 0; 5 < m < 49 and 100 < i < 1000, where m is the total number of particles in 
the grid on the average and / is the total number of occurrences of the grid. 



namely [—90°, 90° ] in both 9 and 4> variables into 7 equal intervals resulting in 49 equal 
sized grids. We then computed the evolution of initially uniformly distributed particles of 
different aspect ratio within the range of rg equal to 0.2 to 2 in steps of 0.2. 

We followed the evolution of the ensemble of particles for 3000 iterations of the stro¬ 
boscopic plot and deleted the first 1000 values to remove transients in the case of constant 
and zero forces and in the case of periodic forces with control and without control. In all 
cases we calculated the number of particles in each grid on every second iteration of the 
Poincare section of the evolution equations resulting in a total of 1000 values. We noted 
the grids in which the total number of particles was greater than or equal to 5 and also 
noted the number of particles in each grid only if the particle occurred in that grid in more 
than 100 iterations in all cases. We denote these values as /, m where /, m denote 
the aspect ratio, total number of occurrences of the grid and total number of particles in 
the grid on the average respectively and prepared tables for these values. A particle of a 
given aspect ratio visiting a given grid and the absence of particles of other aspect ratios 
visiting the same grid is more important than the number of visitations of a given grid and 
the number of particles visiting a grid from the point of view of particle separation. Since 
our earlier computations indicated the greatest sensitivity of the results to the aspect ratio 
near k 2 = 10, we selected k 2 equal to 12.0 for comparing the effects of periodic forces 
with and without control and as a magnitude of the constant force. 













































Figure 7. Stroboscopic plot with control showing a period-5 solution for it 2 = 12; kt 
initial conditions 0 = <p = 45°; co = J;re = 1.6. 





Figure 8. Stroboscopic plot with different controls showing two different period-2 orbits 
obtained for ki = 12; initial conditions 6 = 4> = 45°; w = 7; = 1.2 with control (a) 


it' = 4.0, (b) it' = 5.0. 
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Figure 9. Stroboscopic plot showing period-2 and period-3 orbits successively obtained 
with control applied at every second and third periods respectively for k 2 = 12; initial 
conditions 0 = (p = 45*^; = 7; = 1.6; = 5.0. 


We noted the evolution of the ensemble of particles at every 2nd iteration. We analysed 
in detail the cases of periodic force with and without control and compared the results with 
the constant force and zero forces. The fact that a chaotic behaviour can be controlled to 
oscillate in a selected period resulted in more efficient separation enabling all initially and 
uniformly distributed particles to be directed to a desired grid by a suitably engineered 
control technique. For example, all the 49 initially and uniformly distributed particles of 
aspect ratio = 1.2 could be brought to the grids (1,4) or (7, 4) by applying k!^ = 4. 
Under the same circumstances it was possible to bring all the 49 particles of aspect ratio 
rg = 1.2 to the grids (1, 3) or (7,1) by applying a control of constant force ^2 = 5. A 
sample of the results we obtained for a periodic force of magnitude ^2 = 12 with control 
forces k 2 = 2, 3 and 5 and for constant and periodic forces of magnitude k 2 = 12.0 
and zero forces are given in tables 1-6. The tables presented in this paper indicate that 
controlling chaotic dynamics is preferable to chaotic behaviour and regular behaviour for 
efficient separation of particles. 

A detailed analysis of all the tables indicates that particles of aspect ratio 0.2 alone can be 
separated from a mixture containing particles of different aspect ratios ranging within 0.2 
to 2.0 by applying a control force of magnitude ^2 = 2,3,4 or 5 applied at every second 

period along with k 2 = 12, since particles of aspect ratio 0.2 alone occur in some grids 
as can be seen from the sample tables. In the case of periodic forces without control and 
constant and zero forces particles of aspect ratio 0.2 alone also occur along the boundary 
of the tables. Kumar et al (1995) suggested that particles of this aspect ratio alone could 
be separated by applying a constant force. However, in the case of a periodic force with 



Table 4. Distribution of evolution of initially uniformly distributed particles of different aspect ratios 
for the case k 2 ~ 12; /c^ = 2; ru = J; 5 < m <49 and 100 < Z < 1000, where m is the total number of 
particles in the grid on the average and / is the total number of occurrences of the grid. 
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control, the occurrences of this particle alone are concentrated among a fewer number of 
grids and visit a given grid a larger number of times. In the case of zero forces the particles of 
this aspect ratio alone are concentrated among a fewer number of grids along with particles 
of higher aspect ratio. Thus, for the separation of the particles of this aspect ratio alone a 
periodic force of magnitude k 2 = 12.0 with control is preferable to any other possibility. 
Similar analysis of the tables shows that it is desirable to apply a periodic force with control 
for separating particles of aspect ratio 0.4. As can be seen from the sample tables particles 
of aspect ratio 0.6 alone and 0.8 alone also can separated individually by applying a periodic 
force with control. In the case of particles of aspect ratio 0.8, they can be brought to one of 
the extreme grids. Thus, a periodic force with control seems to be preferable to any other 
case considered for effective separation of particles in the case of aspect ratios rg < 1. 

Kumar et al (1995) demonstrated that periodic forces are preferable to constant forces 
and zero forces for separating particles of aspect ratio rg > 1. In their analysis they 
observed that the occurrence of particles of aspect ratio rg > 1 are spread among a large 
number of grids. However, for particles of aspect ratio greater than or equal to 1.6 individual 
separation was not possible in their analysis. Further study of the tables prepared in this 
work indicate that a periodic force with control is once again preferable in such cases. In 
the case of particles of aspect ratio rg > 1.0 all the initially uniformly distributed particles 
of the same aspect ratio could be brought to one grid except for rg = 1.8. Hence particles of 
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aspect ratio within the range of ranging from 1.2 to 2.0 except for = 1.8 can be easily 
separated as can be seen from the tables. For particles of aspect ratio equal to 1.8 individual 
separation may not be possible in the cases considered in this work since particles of this 
aspect ratio appear in combination with particles of lower aspect ratio equal to 0.2 in some 
cases considered with control. Hence if particles of aspect ratio 0.2 have been separated 
out as explained earlier, particles of aspect ratio equal to 1.8 can be separated from the 
mixture. One advantage of periodic forces with control is that all particles of aspect ratio 
1.8 can be brought to one grid in combination with the particles of aspect ratio 0.2 as 
shown in table 6. 

In conclusion, it has been generally noted that control of chaotic behaviour gives better 
separation than chaotic and regular behaviour. One of the main features of the method 
suggested in this paper is that all particles of the same aspect ratio can be concentrated in a 
previously desired grid. A detailed analysis of the problem can suggest suitable designs for 
separation of particles by aspect ratio to get a well-characterized suspension of particles. 
A possible design for this separation of particles from a mixture of particles of different 
aspect ratios based on the differences in the orientation of the particle may consist of a base 
plate having grooves along different orientations so that when the particles are oriented in 
such directions, they settle in a particular groove and can be separated out at every integral 
multiple of the fundamental period. 
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Thus, even in the rather simple application considered in this work, control of chaos 
leads to greater efficiency. This suggests that the possibility of chaos control should be 
important in many of the applications mentioned in the introduction. The novel control 
of chaos technique suggested in this paper has been demonstrated to be effective even in 
the relatively complex problem considered here. An additional feature of the control of 
chaos technique suggested in this paper is that the control is effected very rapidly and the 
behaviour of the concerned system can be switched from one desired period to another 
desired period very quickly. This suggests that this control of chaos technique may be 
applied to other chaotic systems very effectively. 
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b equatorial radius; [L] 

i unit vector in the x direction; [-] 

J 2iz{re + Ve'^) [-] 

k periodic force field; [MLT“^] 

k' constant force field; [MLT~^] 

intrinsic properties of the particle, which depend upon 
the geometric configuration of its wetted surface; [L^] 

^K, T intrinsic properties of the particle, which depend upon 

the geometric configuration of its wetted surface; [L^] 

I total number of occurrences of the grid; [-] 

L torque about the centre of the particle exerted 

by the fluid on the particle; [ML^T~^] 

m total number of particles in the grid on the average; [-] 

n period of controlling; H 

re particle aspect ratio given by = a/b', [-] 

5 rate of deformation tensor given by S = ^[VV + VV^]; [T~'] 

u* unit vector along the axis of symmetry which determines 

the orientation of the spheroids with components Ui*, U2*, U3*; [-] 

U translational velocity; [LT~^] 

U —V translational slip velocity; [LT“^] 

y y coordinate; [L] 

y shear rate; [T“^] 

jjiQ viscosity of the fluid; [ML~ ^ T~ ^ ] 

6, <p polar angle and the azimuthal angle corresponding to 

a given vector u respectively; [-] 

6 time derivative of 6>; [-] 

0 time derivative of 0; [-] 

0 ) frequency of the driving force; [T“^] 

O angular velocity of the particle; [T“^] 

fi — CD rotational slip velocity; [T“^] 

V del operator. [L~^] 
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Abstract. This paper investigates the non-terminating (NTB) behaviour of 
discrete-event processes (DEPs) modelled by finite state Biichi automata. 
Specifically, we develop synthesis criteria and procedures for state-feedback 
supervisors that ensure infinitely often (i.o.) returns to a given good state sub¬ 
set of the DEP in the closed loop. This problem is referred to as the stabilization 
problem for discrete event systems (DESs) elsewhere. 

The approach of the present paper is distinguished by the fact that the stabi¬ 
lization is achieved without major compromises on the richness of behaviour, 
or manouverability. This flexibility comes as a consequence of using two super¬ 
visors in parallel and switching them on-line to achieve the desired behaviour 
characteristics in a non-blocking manner. Supervision is arrived at by decom¬ 
posing the state space into controllable and uncontrollable subspaces relative 
to the desired target space, much in a manner identical to control theory for 
systems modelled by differential equations. 

Keywords. Regenerative supervision; non-terminating behaviour; discrete 
event processes. 


1. Introduction 

Supervisory control of discrete-event systems (DESs) modelled as finite-state machines 
(FSMs) is an important control research area presently. The basic model and control theory, 
in an FSM-Formal Language set-up, are due to Ramadge and Wonham (Ramadge & 
Wonham 1987; Wonham & Ramadge 1987). Initially, the supervisory theory was developed 
only for the transient behaviour control of a DES, which behaviour arises, naturally, from 
simple FSM models and the regular finite-string behaviour sets of such models (Eilenberg 
1974; Hopcroft & Ullman 1979). Later the theory was extended by Ramadge (1989) and 
others (Kumar et al 1992) to control the steady-state DES behaviour by using Biichii 
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extensions of FSMs, or the so-called Biichi automata (Eilenberg 1974), which model non- 
terminating (or sequential) response of the FSMs. 

An important problem that arises in the study of sequential DES behaviour is that of 
stability. Ozveren & Willsky (1991) define stability as infinitely often (i.o.) guaranteed 
returns, by means of a state feedback, to a set of a priori selected good states of the DES, 
in a non-Biichii setting. This paper addresses a similar problem in a Biichii setting which, 
as mentioned above, arises as the natural extension of FSMs (Eilenberg 1974) to model 
sequential behaviour. Furtheimore, the control strategies developed here ensure flexibility 
in the sense that the rate of convergence to the good set can be varied on-line, by adjusting 
a single free parameter, so that depending upon other (possibly unmodelled) conditions 
one can trade-off richness of behaviour (manouverability) for quick stabilization. This 
is done by having two logically connected parallel controllers which are switched based 
on the requirement at any given instant during the operation or evolution of the process. 
Moreover, the controllers we develop satisfy the required mimmal-restrictiveness and 
non-blocking conditions (Ramadge & Wonham 1987), so that the closed loop behaviour 
is maximally-permissible (Wonham & Ramadge 1987) and satisfies liveness conditions 
(Ozveren & Willsky 1987). Another important feature is that finite and sequential - that 
is transient and steady-state - specifications are simultaneously handled, thus simplifying 
the overall control process. In addition, the stringent requirement of topological closure of 
the specification behaviour (Ramadge 1989) is avoided here. 

There are other notions of stability as well. For example. Brave & Heyman (1989) 
consider the good set as an attractive invariant domain in the closed loop under state 
feedback. The control problem then is to design a feedback so that the DES enters the given 
set in a finite number of steps and does not exit it later. This could, obviously, constrain 
the closed loop behaviour seriously since it compromises the richness of behaviour. Other 
notions of stability that involve the sublanguages of the behaviour of a DES (Kumar et al) 
rather than state subsets of the underlying FSM also exist. We, shall, however, deal only 
with the state-set defined case, and call onr flexible stabilizers regenerative supervisors, 
to highlight the fact that the controller ensures, in the closed loop, i.o. visitations to the a 
priori chosen good state subset. 

This paper is organised as follows: § 2 develops the necessary formal background for 
DEPs and their supervisors in terms of both state dynamics and language behaviour. The 
main section of this paper is § 3, which lays down the synthesis theory for design of the 
required flexible supervisors. An important part of this section is the proof of the necessity 
condition of theorem 1, which is separately given as appendix A due to importance beyond 
the needs of the theorem. Section 4 concludes the paper by highlighting the advantages of 
our procedure with the help of a near generic example given in appendix B. 


2. Mathematical preliminaries 

2.1 Languages 

Let S be a finite set and S* the set of all finite strings over S. E* is a free monoid (Eilenberg 
1974) under catenation with the empty string e being the monoid identity. For 5 e S'*', h i 
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(proper prefix) of v, denoted hy u :< v (u < v), if for some co e E*, v = uco i\co\ > 0). 
Clearly the relation is a partial order on E*. 

A subset L c S* is called a Language over S, and the set L = {u\uv e L for some v € 
E*} is called the prefix closure of L (Eilenberg 1974). L itself is prefix closed if L = L. 
Next we define cu-Languages (Eilenberg 1974; Ramadge 1989) which constitute infinte or 
sequential extension of strings. 

Let N denote the natural number set and let S" denote the set of sequences over E, i.e., 

= {e\e\N —> S}. (1) 

For e and j e N, let e(j) denotes the yth element of e, and e^ = e(l)e(2 )... e(j) 
denote the finite string (i.e., e E*) consisting of first y elements of In general, m e Y* 
will be called a prefix of e € E'" (denoted, once again, by u :< e) if u = for some 
7 e A. 5 c E" is called an <y-Language over E, and the language pr(B) g E* defined 
as; 

pr(B) = {€}U{ej\ e e BJ >1}, (2) 

is called the prefix of B. 

For a monotonic sequence ui < U 2 < m-- of elements over E*, 3 a unique element 
e G E" such that e^ = uk, for j = \uk\,k G N (Eilenberg 1974). Such an e is called 
the limit of the sequence {uk}. The adherence or limit of a prefix-closed L c E* is the 
du-Language defined as; 

L°° = {e\e G E" A e^ G L, V; G A). (3) 

If L is not prefix closed it suffices that e^ G L for infinitely many j € N (Eilenberg 
1974; Ramadge 1989). 

A crucial metric topology on E" is induced by the metric p, defined as under (Eilenberg 
1974); 

P(^l) ^2) = if and e" 7^ e^', and, (4) 

p(ei,e2)=0ife^ =e2. (5) 

It is easy to show that in this topology Y^ is compact (Eilenberg 1974). Furthermore, 
this metric can be extended in a straightforward manner to E* U E*", though we shall not 
need that extension here in the sequel. The topological closure B (B union its limit points) 
of an d6)-Language B equals the limit of its prefix language, i.e., B = [pr(J5)]°°, as can be 
readily verified (Eilenberg 1974; Ramadge 1989). For 5 c 5 c we say B is closed 
relative to S if B D S = B (Ramadge 1989). 

2.2 Discrete-event systems 

A discrete-event system (Ramadge & Wonham 1987) is modelled as a finite state au¬ 
tomaton (Eilenberg 1974) G = (Q,Y, 8, qo) where Q is the (finite) state set, E the 
finite set of events (also called the process alphabet), 8 \ Q x Y —^ Q gives the 
state transition (partial) function, which can be naturally extended to E* by defining 


8(q,€) = q, and 8{q, sa) = 8(8(q, s'), cr) Vs € S*, whenever all intermediate transi¬ 
tions are defined; else it will be undefined. Whenever 8(q,s) is defined for some s e S* 
we shall denote it by (5(^, s)!. 

We identify the following subsets of S and Q: 

VqeQ- : <5(^,u-)!}, (6) 

e (2; VX C Q; {^t E S ; 8(q, a) £ X), (7) 

VX c (2; TX:=[q eQ : Bct E E s.t. 8(q, cr) E X). (8) 


Eg defines the active event set at q, i.e., set of events executable at rX is the source 
set of X c Q, that is the set of states which enter X in a single step and E^^x ^ E^^ is the 
active event subset at q which leads into X. For all state subsets X c Q, we shall denote 
by a(X) the cardinality of X and by X^, the complement of X, In particular, a(Q) < oo 
and 

For the DES G the finite (transient) behaviour will be given by the prefix closed lan¬ 
guage L(G) = [w\w E S* and 8{qQ, in)!}. Sometimes, in order to emphasise that this 
constitutes the event labels of the paths starting at qo of the underlying automaton G, we 
may, instead, call this behaviour L(G)^q. We shall also mean by the natural sequential (or, 
non-terminating) behaviour (NTB) of G the limit [L(G)]°^ which will be assumed to be 
non-empty throughout. In the simplest form the Biichii extension of G can be defined as 
follows (Eilenberg 1974; Ramadge 1989). Let Qm c g be a given state subset. For each 
e E [L(G)]°°, 3 a unique Se'. N —^ Q, such that Sdj) ~ 8(qQ, e^), j E N. The state 
sequence Se is called the state trajectory for the behaviour sequence e. The sequence e 
and the state trajectory Se are said to be admissible, relative to Qm (or, sim^ply admissible) 
if Se visits the state subset <2m i-o., i.e., if Se(j) e Qm for infinitely many j E N. The 
Biichii set Qm then defines the acceptable sequential behaviour of the DES G as follows: 

[L(G)]°°(2m = 1 ^ [f-(G)]°° ) and Sg are admissible). (9) 

In principle, one can use multiple Biichii (Eilenberg 1974) sets to generate very complex 
behaviour. In the sequel, however, only one set will be used. Furthermore, it is obvious 
that [L(G)]°° itself can be taken as the Biichii behaviour of G for Qm — Q. 

The control mechanism on DES G = (g, E, <5, ^o) is defined as follows (Ramadge & 
Wonham 1989): on E a partition E = E^ U Ej, is defined, where the events in E^; are 
called the controllable events and those in E„ the uncontrollable events. Let E := (0, 1}^"' 
be the set of binary assignments on Ec, i.e., the discrete topology on E^. Then each 
y E r, (y : Ec —^ {0,1}) is a control pattern, in the sense that we say <t e Ec is 
enabled if y(cr) = 1 and is disabled otherwise. Clearly y e P can be extended to whole 
of E by defining y(cr) = 1, Vcr e E^. Given a DES G = (Q, S, 5, q^) , the controlled 
DES Gc := (C» El, (5c, ^o) is defined as follows: 


5c: r X (2 X E* Q s.t., 

if (g,or)! andy(or) = 1, 
I undefined, otherwise 


( 10 ) 


Throughout this paper we shall assume that the E partition is given and the DES admits 
the above mentioned control mechanism. We shall, henceforth, drop the subscript c and 
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assume G — (Q, I], 5, q^) to be a controlled DES, or a controlled discrete-event process 
(CDEP). 

2.3 Supervisors for CDEPs 

There are two main supervisory control schemes for CDEPs: (a) dynamic (output) feed¬ 
back, and (b) state-feedback supervisors. 


2.3a Dynamic supervisors: A dynamic supervisor (Ramadge & Wonham 1987) for a 
CDEP G is a function M : L(G) —P, which specifies for each string s e L(G) the 
subset II5 c E of enabled events (E^ := {a : M{s)o = 1}). The closed-loop behaviour 
of G supervised by M, denoted by L(G/M), is given by: 


€eL(G/M), and 


o)cy G L{G/M) -4^ 


(0 

(ii) 

(Hi) 


CO E L(G/M), 
M(co)a = 1, and 
(oa E L(G). 


( 11 ) 


The non-terminating (or, steady-state) behaviour of G supervised by M is, obviously, 
then given by [L(G/M)]°°. 


2.3b State-feedback supervisors: More important, from the point of view of applica¬ 
bility and design simplicity, is the state-feedback supervisor (Kumar et al 1993a) which, 
for a given CDEP G is defined as a function R : Q —^ E. Such a supervisor speci¬ 
fies, for each state of CDEP, the set of enabled events. The closed-loop system, where G 
is supervised by R, can be viewed as another state machine (G/R) := (Q, E, 5/?, ^0), 
where 5/?(<7, cr) = 8(q, a) if R(q)<J = 1, and, 8R(q, a) is undefined if R(q)a = 0. The 
non-terminating behaviour of the closed-loop is, once again, given by [L(G//?)]'^. 

An important property that any supervisor should satisfy is that every response sequence 
begun (in the closed-loop situation) must be allowed to reach a steady state; in other 
words, every finite string generated by the closed loop must lead to a non-terminating 
behaviour sequence. This is called the property of non-blocking or the liveness property 
of the system which implies that the system does not stall along a response trajectory. 
Formally, a supervisor M (dynamic or static) for a CDEP G is said to be non-blocking if: 

Ms E L(G/M)\ 3e G [L(G/M)]°° s.t. s — e^ for some j e N. (12) 

Given a CDEP G, and an X c g, a supervisor M is said to be X-invariant if: 

Vs G L(GtM), S(qo, s) E X. (13) 

An X-invariant (or, equivalently, a state-invariant) supervisor will, therefore, restrict the 



[L(G/M)]~ 7 ^cJ), and 

[L(G/M)r c[L(G)r2,„- (14) 

The supervisor M for G is such that the NTB is non-empty and the set Q,n is visited i.o. 
Note that this does not imply that the only states visited i.o. are the Qm states. Also note 
that this is identical, in principle, to the stability notion of Ozveren & Willsky (1991), as 
introduced earlier, although the setting here is Biichii in contrast, wherein the sequential 
behaviour is defined in terms of the admissibility conditions as given in (9). 

The problem addressed to in the sequel is the following; 

Given a CDEP G = (Q, S, 5, qo) and the partition S = U Sc, and a state subset 

Qm ^ Q\ 

• Determine if there exists a supervisor M for G that is both non-blocking and Q,n- 
regenerative. (Henceforth, we shall simply call these supervisors regenerative once the 
set Q,n is understood as identified). 

• If Ma is the non-empty family of non-blocking regenerative supervisors then find M e 
Ma such that: 

[L(G/5)]°^ c [L(G/M)]°^ G (15) 

i.e., find the minimally restrictive of all the non-blocking supervisors for G. 

We shall call these supervisors simply regenerative whenever the Biichii set Qm is 
specified or obvious. 

3. Main results 

In this section we develop a controller synthesis for solving the above mentioned stabiliza¬ 
tion problem. The synthesis procedure is based on construction of a series of state subsets, 
which are recursively reachable and within the union of which the state invariance of a 
supervisor can be easily achieved. 

3.1 Control subsets 

Given G = (Q, E, 5, < 70 ) , S = Eu U Ef, and ^ Q, define the state subset 

sequences Qi and A/ as follows: 

Qq = Xo = Qm, (16) 

Xi^\ = {qeiQi^CtQi): c E^}, (17) 

Q/+]~( 2 /UX,-+i. (18) 

Note that Xi , for example, is that subset of Qq which can be driven to Qo in a single tran¬ 
sition. This is because there is a transition to Qo, and all transitions to Qq are controllable 



(and, therefore, can be disabled). Since {} is an ascending chain we define. First Control 
Subset (FCS) Q as its limit, given by: 

(2= sup <2/ = 

/ i 

The FCS Q is, therefore, the largest state subset which reaches Qm, possibly uncontrol¬ 
lably, but whose dynamics outside of itself, being controllable, can be disabled. To control 
escape dynamics from Qm itself, we define. First Regenerative Set (FRS) (2,„ as: 

Qrn = {q^ iQm H TQ) : c E,}. (20) 

The state set Q,^ is clearly the Biichii subset which can be made invariant in the closed 
loop by disabling the exit transitions. 

Some interesting properties of the sets X/, Q and which essentially identify con¬ 
trollable transitions at various levels, are listed below (see Valavi for more properties and 
proofs): 


3k < a(0^^) s.t. Xi = 4> Wi > k, and, 

Mi < k, Xi ^ <J>, and 

k 

Qk=Q = 

(21) 

i=0 

Vi >l,Vg€ Xi, IJ C Sc. 

(22) 

j>i 

Mq E Q — Qq, ^q^Qc ^ 

(23) 

/-I 

Vi > 1,V? € Xi, IJ 

(24) 

j=0 

Intuitively, these properties all stem from the structure of the Xi 

chain, which is such 


that all dynamics from the lower to higher index can be disabled. Therefore, as the first 
approximation, the supremum of the chain contains all the states that can be driven to Qm 
either uncontrollably, or by enabling all controllable events along the path. Also note that 
property 21 is valid for a finite Q, whereas others hold independently of its cardinality. 
We shall, however, consider our CDEPs to be finite state only, as earlier mentioned, for 
practical reasons, even if some of the results hold in general. 

Next, we develop the regenerative supervisors for our problem as defined in §2.3. This 
will be achieved by synthesis of a Q-invariant supervisor, so that all potential paths that 
escape Qm indefinitely are disabled. As will be clear from the sequel, it is always possible 
to do so under the fairly mild condition that be in Q , i.e., in one of the Xi ’s, which are 
state sets that can always be driven to Qm in a finite number of steps. Initially, to to keep 
the synthesis simple and tractable, we shall assume that the FRS Qm equals Qm • Later 
this condition will be relaxed. 

Note that V( 3 ’ € Qm — Qm‘^ either E^ ^ = OorE^ is uncontrollable. For regenerative 
behaviour, it should be possible to disable transitions to Q^, and for non-blocking condition 



there should be a transition to Q. Therefore, Qm — is the set of states that violates one 
or both of these conditions. After the supervisory syntheses given below this problem will 
be addressed to, so that the final synthesis is free of this apparently severe constraint that 
Qm equal Qm- 


Theorem 1. {(Jhe R-siipei-visor)] Given the CDEP G = {Q, E, 5, qo), with S = U 
Sf: Qm Q Q> the PCS Q and the FRS Qm, assume that Qm = Qm. Then there exists a 
non-blocking and regenerative supervisor Rfor G if and only ifqo e Q. 


Proof. Sufficiency: Let qo ^ Q. Then from the above properties we know that Q = 
ULo h < q:(<2^„) < a(Q) and Xi ^ <i>, 0 < i < k. 

Define a state feedback supervisor R as follows: 

Wq e Qo, R(q)cT=0\/a € E^^^, (25) 

[c E,.; (20)], 

Vi>l,^qex,, i?(9)o' = 0Va€(U)>iS,,Xj)U(S^^,), (26) 

[c Ec; properties 21, 23] 

R(q)a = I, otherwise. (27) 


R, as is obvious, is a supervisor that drives the system fastest to Qm, the desired state set for 
regenerative operation, by preventing escape from Q as well as any backward movement 
(towards increasing index on Qi) within Q. Therefore, once the system is in Q, it will 
move into Qm in at most k steps. 

The closed loop system G/R = (Q,I], Sg, qo), is then defined by: 


^Dr /7 n-'i - if i?(<?)o-= 1, and 

oRKq, jundefined, if R(q)a = 0 


(28) 


^ To show that/? is note that 6 Q, 5^(^, a)! (isdefined)=^> SR{q,a) £ 

Q, because transitions from Q to its complement are disabled by the supervisor R ((25), 
(26)). Therefore, it follows that: 

qo e Q ==^ Sg(qo, s) G Q, Vs € L(G/R). (29) 

In other words, /? is a ^-invariant supervisor. Furthermore, (20), (24), (27) and the fact 
that Q^ QoG (Uf=i ^i) together imply that: 

Vq gQ\ {a G E^^^ : R(q)cT = 1) ^ (I>. (30) 

Next, from (29) and (30), it follows that € L(G/R), Bcr € E suchthat^a e L(G/R). 
Therefore, V^o 6 L(G/R), 3 a sequence 5o ^ ^ ^2 • • * in L(G/R), so that R is sl 

non-blocking supervisor that keeps the closed-loop live. 

To show that R is regenerative, it suffices to show that a transition to Xo(= Qm) is 
guaranteed in at most k steps. To this end assume that sco g L(G/R), with 8R(qo, s) = 
q G Xr, for some i < r < k (this is guaranteed by <^0 ^ Q ^d the supervisor design 
(25)-(27), and |a>| > k. Let crjt-icryt _2 • • • cro be the k-long prefix of o) (note the reverse 
index on co). Clearly qj = 8R(qj_i, Ok-j) e Xrj, where qo = q and r > ri > r 2 • • •. 



Since the index r/ decreases at least as fast as the index on cr, we have qi e Xq for some 
i < k. Furthermore, if ^ e Xq, 8R{q, cr^_i) e Xq\J X,-, for some 1 < r < /:, the same 
process holds. Thus, in the closed loop a transition to Qm{= Xq) is always guaranteed (in 
< k steps) and, therefore, the supervisor is regenerative. 

Necessity. See appendix A. □ 

Remark 1. It is obvious from appendix A that regenerative and non-blocking supervisors 
are ^-invariant, since for every sub-CDEP that evolves from a reachable state q e Q, 
we have that q & Qi (for some i) Q Q. This implies that only the states from Q can be 
regeneratively driven to Qm along non-blocking trajectories. 

3.1a Minimally restrictive supervisors: The /^-supervisor, though non-blocking and 
regenerative is, very restrictive in terms of the richness of the closed-loop behaviour. In 
fact, as is obvious, the /^-action is to select control patterns to drive the closed-loop to Qm 
as fast as possible. This constrains the behaviour profile of the system. To overcome this 
problem we use another supervisor in parallel with the /^-supervisor, so that the closed-loop 
is both rich in behaviour {minimally restrictive) as well as regenerative. 

Define a state-feedback supervisor 5 : Q —> F for the CDEP G as: 

'iq € < 2 ; 5(g)or = 0 Vcr € [C S,], (31) 

S{q)G = \ otherwise. (32) 

li qQ e Q the supervisor 5 is 2-invariant, since e Q the transitions that lead to a 
state in are disabled. Note that the supervisor S, while ensuring that the system remain 
inside Q , does not drive it hard to Qm, the way R does. We now have the following 
straightforward result. 

Theorem 2. Let G be a CDEP, as in theorem 1, with Qm = Qm> qQ ^ Q. Let S be 
the state feedback supervisor for G as defined above. Then, if M is any non-blocking and 
regenerative supervisor for G, we have that: 

L(G/M) c L{G/S). (33) 

In other words, S is the minimally restrictive ^-invariant supervisor, richer in behaviour 
than the non-blocking and regenerative supervisor subclass. 

Proof (Straightforward) We know that every non-blocking and regenerative supervisor 
is <2-invariant under the given conditions (see remark 1). Therefore, since qo Q Q and 
the states in Q^ cannot be driven to Qm in the desired manner, every regenerative and 
non-blocking supervisor must disable transitions to Q‘^. 

Therefore, .v 6 L{GIM) 8{qo,a)) e Q, Va> e pr{s). Also since the only transi¬ 
tions that S disables are those to <2^» we have that 8{qQ,a)) e Q ==> a) G L{G/S), \fo) e 
pr{s) s G L{G/S). The empty string e in pr{s) does not affect the argument since 
qo G Q. This completes the proof. □ 


3.2 Flexible supervisors 


From the preceeding discussion it is obvious that the /^-supervisor capability of fastest 
motion to the desired state set constrains its behaviour, whereas, in contradistinction, 
the 5-supervisor’s minimally restrictive behaviour limits its frquency of arrival into the 
desired state set. In practice, a compromise is desirable. In the rest of this section we shall 
address this problem and provide a solution based on continual switching between R and 
5 supervisors that, dynamically, takes into account the last arrival into Qm. 

First we define a delay function g : L(G) —^ N U {0} (N is the set of natural numbers) 
as: 


5(e) = 0, 


g{SG) = 


g{s) -f 1 
0 


if ^{qQ,sa) 6 (2^, and 
if 8{qo,sa) e Qm- 


(34) 


The delay function is aperiodically cyclic and in every cycle counts the number of 
events executed outside of Qm, so that a large value in any cycle reflects large behaviour 
executed without returning to Qm- Next, given a. k g N, we. define a k-jump function 
/k:/VU{0}-^{0, Ijas: 


fk(i) = 


0 , 

1 , 


if i < k, and 
if i > k. 


(35) 


g fj^ 

The composite sequence L(G) —^ N U {0} —^ {0,1}, represented by the function 
composition /k o g, therefore, becomes high every time the CDEP executes more than k 
events in a row outside Qm- The composite function /k o g will henceforth be called the 
system switching function of order k. 

The switchingJun.ction is used for control strategy as follows: Run the CDEP G under the 
supervisor 5 as long as the switching function is low (= 0) and once the switching function 
becomes high (= 1) the supervisor R takes over. The parameter k can be dynamically 
varied to satisfy the requirement on the rate of visitations to Qm - Note that k = 0 implies 
running the system only under R, since the fastest convergence to Qm is desired and as 
—> oo only 5 will be used since the constraint now is the richness of behaviour rather 
than fast settling to Qm- The composite overall dynamic supervisor (see 2.3a) T modelling 
the abovementioned control strategy can be defined as under. 


Definition. Given a CDEP G = {Q, S, 5, qo), E = U and the associated state 
sets Qm, Q, and Qm such that Qm = Qm and qo G Q. Let A: € N U {0} be given and the 
supervisors R and S be as defined earlier. Further let the switching function /k o g be as 
defined above. 

Define the dynamic supervisor T : L(G) —^ F as: 

T(s)g = {[/k o g(j')] A [R(<5(^0. ‘y))or]} V {-’[/k o g(s)] A [5(<5(^o, ‘y))o']}. 

(36) 

From the preceeding discussion it follows that for finite k the supervisor T is both 
non-blocking and regenerative with the richness of behaviour increasing as k becomes 



ly specification changes or needs of the environment. 


The computational effort involved in T synthesis is mainly in the recursive computation 
f the sets X/ , so that a loose (conservative) upper bound on the complexity is Q){n^p), 
^here n = ci{Q) and p is the number of edges of G. 

.3 Weak supervisoiy conditions 

lS explained in § 3.1, for regenerative behaviour it should be possible to disable transitions 
) and for non-blocking behaviour transitions to Q must be ensured. The states in 
— Qm precisely those Biichii states that violate one or both of these conditions, 
herefore, for desired behaviour it must be ensured that these states are kept out of the 
ossible closed-loop paths with minimal reduction in the desired behaviour set. 

This can be done as follows; Let G be given such that ^ Qm- We proceed to build 


family of control subsets as follows. 

Qmo — Qm\ Qo ~ Qmo ~ Qm - Define! 

Qmi = Qmi-\ ’ 

Qi = set obtained from PCS construction with 2m,.as Q^, (38) 

set obtained from eqn. 20 with Qim^^ Qm^tid Qi as Q (39) 

lote that the chains Qmi , Qi and (2,n/ descending. Next define: 

Qs = mfQi = 

I . 

I 

2m, = inf e„, = HQmi, (41) 

I . 

I 

Qm, = inf <2m; = r\Q>nr 

i 


The role played by Qs and is the same as that of Q and Q„i, respectively, in the ear- 

er case when Q was assumed equal to . In other words, transitions to can always 
e guaranteed from Q^, whereas from outside of Qs this cannot be achieved. Therefore, 
/e shall call Qs the Infimal control subset (ICS) and the Infimal regenerative subset 
IRS). It is easily verified that — Qm, • We call the set Qm, the Infimal marked subset 
IMS). It is obvious that if the algorithm of theorem 1 is ran using Qm as Qm, Qs Q 
nd Qm^ as Qm, then the following generalization of that theorem holds. 

theorem 3. (Generalization) Given a CDEP G = (Q, S, (5, go), with E = U Tic, 
nd the Biichii state subset Qm ^ Q, Qs and Qm, t/ie associated ICS and IRS, 
espectively. Then G is a regenerative CDEP (Le., there exists a regenerative supervisor 
or G) if and only ifqo € Qs- 


The example given in appendix B highlights the advantages of the synthesis given here 
relative to the well-known supervisory synthesis of Ramadge (1989). The principal feature 
in favour of synthesis in by Ramadge (1989) is that the NTB specifications are given in 
terms of subsets of [L(G)]°° , which is more general than specifying i.o. visitations to a 
given Biichii state subset. However, the approach requires conversion of the NTB to finite 
behaviour using the prefix operator pr(.) and then synthesising the supervisor that meets 
the prefix specifications. The limitations of this method are as below. 

• A supervisor having a finite state structre is possible only if pr{B) is regular (Ushio 
1993). 

• The synthesis requires that B be closed relative to [L (G)]°°, which may be hard to obtain 
in practice (see appendix B). In case B does not satisfy this constraint, the closed-loop 
can only be guaranteed to stay within 5, the topological closure of B relative to [L(G)]°^, 
and not in B which is the NTB specification (see appendix B). 

The acceptable specifications used here, though relatively less general, are meaningful 
and adequate for a large CDEP class. Main advantages of our method, for example, are 
the following. 

• The Biichii set specification is veiy meaningful in teims of closed -loop stability analysis 

(Ozveren & Willsky 1991). In fact, stability specifications in terms of Biichii sets is a very 
natural and logical extension of CDEPs based on finite state machine models (Eilenberg 
1974). The synthesis then provides stabilization of CDEPs. 

• The method does not require construction of a grammar for acceptable NTB, thus re¬ 
ducing the computational complexity. 

• It is easy to deal with, possibly conflicting, finite and non-terminating behaviour. This 
is due to the supervisor structure obtained and the operational flexibility achieved due 
to our switching junction. 

The supervisory synthesis approach as outlined above, therefore, presents a very effi¬ 
cient mechanism of ensuring flexible stabilization and satisfaction of other closed-loop 
specifications of CDEPs. 


The authors are deeply grateful to one of the anonymous reviewers for pointing out a 
serious error in the proof of theorem 1 as well as other minor errors. They are also thankful 
to other reviewers for many helpful suggestions. 

Appendix A 

Completion of the proof of theorem 1 requires some auxiliary definitions and lemmas. 
We begin by introducing the notion of a restricted supervisor. Let G = (Q, S, 5, ^o) t>e 
a CDEP. Suppose <5(^0, s) = q e Qfor some s € E*. Then G = (Q, E, q, S), where 



Biichii set for the 5 mZ?-CDEP will be given by Qm = m • 

Definition. Let G be a CDEP and G a ^wib-CDEP of G. Suppose T ; Q —^ F is a 
supervisor for G. Then f : Q —F is said to be the restriction of F to (2 if: 

f{q)a = T(q)a 'iq € Vcr e S. (Al) 

It follows directly from the above definition and the notions of non-blocking and regener¬ 
ative supervision introduced earlier that if F is a non-blocking and regenerative supervisor 
for G, so is F for G. Next we explore the structure of X, , Qi and Q ((16)-(18)) sets of G 
when the Biichii set is restrictively altered. 

Lemma AL Let G be a CDEP with Qm its associated Biichii set. Let c and let 
Qi, Xi and Q be the control subsets of G relative to Biichii set 2m- Then Qi c g,-, VF 

Proof (By induction on i ). At i = 0 , we have that Xq = Qq = Qy^. Now 2o = Qm =>• 
So ^ Go since 2m ^ Qm- 
Assume that the hypothesis holds for i = j. 

Now by definition (see (17)), we have that q e Xj+] ==^ ^ s (t Qj fi Qj) and S xc C 

Sc. Therefore <7 e r 2i H (2j U (Qj - Qj)) = (r Qj n QpU (r Qj n (Qj - Qj)) because 
Qj c Qj and 2 ; = 2 / U (Qj -Qj). Also note that S xc c Sc c Ec since 

We, therefore, have that Xj+i = Xj+i U € rQj G (Qj - Qj) : £ Sd, 

where the second set is obviously in 2; - Therefore, Xy+i C Xj+i U Qj = Qj+i, and, 
consequently, 2;+i = Qj G Xj+i c Qj u Xj^i ~ 2;+i- Hence Qj £ Qj, V 7 6 N. O 

Remark Al. It is obvious from the above proof that lemma Al holds if G is replaced by a 
5 mZ 7-CDEP and the Biichii subset is compatibly chosen, as discussed earlier. 

Next the proof of the necessity of theorem 1 follows: 

Proof (Necessity of theorem 1). LetG = (Q, S, 8, ^ 0 ) be the given CDEP with Qm(i=- *1^) 
as its associated Biichii state subset. Let F be a non-blocking and regenerative supervisor 
forG. 

We shall prove the necessity of the theorem 1 by induction on 01(2)• 

For Q!(2) = 1, since ^0 = Qm, the result is obvious. Assume that the result holds for 
ci(Q) ~ n. Next, let q;(2) = n -b 1. Suppose 3or 6 S such that ^((^o, or) = q^. Then if 
O' E Sm, we shall have qo € Qm, otherwise no regenerative supervisor can exist as the 
CDEP could indefinitely stay in q^. If this were indeed the case then qo ^ Q and there is 
nothing further to prove. Therefore, without loss of any generality, we shall assume that 
in such a case cr e Ec. 

Next let • • •, < n be the states, other than qo, which are 1-step accessible 

from qo, i.e., Vi < r, qi = S(qQ, crj) for some a,- E E. Consider the sub-CDEPs G/, 1 < 



i S f, generated from the states qi, i.e., Gi = (Qi, E, qi, 8), where Qi c g is Ae state 
set accessible from and Qi,n = Qm^ Qi is the associated Biichii set. Since q!((2/) < n, 
we have, from lemma A1 and remark Al, that, V/, j; / < r and j < ri, Qij ^ 
the jth control subset corresponding to the ith ^m&-CDEP is contained in the jth control 
subset of G. Here ri denotes the length of the control subset chain upto its saturation (see 
( 21 )). 

Furthermore, since the restricted supervisors 7)- acting on G/ are non-blocking and 

regenerative, and since oi(Qi) < n , it follows by induction hypothesis that qi e Qii= 
Qin, by (19)), Vi. Next let r = max/ r/. Then it follows that {q\,q 2 , ■ ■ • ,qr] ^ Qr- 
Clearly then, qo ^ rQf^ D Qj and (since the only such transition can be to 

<?o itself and that is controllable, by assumption); so that q^ 6 Qp+i Q Q. This completes 
the proof. □ 

5. A simple synthesis example 

Consider the CDEP given below in the figure. Here Q ~ [qo, q\, qm}, E = {cri, < 72 , 0 - 3 , 0 - 4 } 
with Sc = { 0 - 2 } (say). Clearly L(G) = o'ia' 2 *(l + crso-;) and [L(G)]°^ = {cricr 2 ", aicr* 
(T 3 cr"}. Let the acceptable (regenerative) non-terminating behaviour (NTB) be infinitely 
often i.o. visitations to the set Qm = [qm] - The corresponding largest desired behaviour 
then is 5 = {cricr^o-scr^} C [L(G)]°^. With this as the desired NTB, the synthesis pro¬ 
cedure of [R] gives a supervisor M, whose finite and NTB behaviour equals L(G/M) = 
cTicrl(1 + and [L(G/M)]°° = {cricr^, aia^cy^cr^]- The closed-loop NTB obtained 

is actually the (topological) closure (see §2.1) of B relative to [L(G)]°° and violates the 
specification of the desired NTB. This limitation of the procedure occurs because B is not 
closed relative to [L(G)]°°. 

If however we choose B = {o-icr 3 or"}, which is closed relative to [L(G)]°°, as the de¬ 
sired NTB, then the algorithm of by Ramadge (1989) results in a supervisor M such that 
L{GJM) = {<y\<y'i<^l] and [L(G/M)]°° = {c 7 ior 3 or 4 "} which satisfies the NTB specifica¬ 
tion. Note, however, that the resulting supervisor turns out to be very restrictive. 

To obtain a less restrictive supervisor, we could choose Bk = {cricr 2 < 73 (j^}, k < 
00 , which is closed w.r.t. [L(G)]'^. With this B the synthesis by Ramadge (1989) 
results in a supervisor less restrictive than the earlier one (larger k improves the per¬ 
formance), but the procedure requires construction of a grammar for pr(Bic) for every 
desired k. 

Using the algorithm of §3, we obtain the following contol subsets as follows. 

^0 = Qm = l^m); Al = {^1}; A2 = {qQ\-Qs — Q'^ Qmg ~ Qm- 


<72 <74 



90 5^1 9m Figure Bl. A CDEP example. 



Observe that qo € Q^. We, therefore, get: 

The supervisor R = R(q)cr z= 0 if q = qi and cr = 02 ; and R{q)<y = 1 otherwise. 

The supervisor S = S{q)<j = I'iq ^ Q and a € S. 

The NTB of the closed-loop under the composite supervisor T is then equal to Bk = 
{cricr 2 (T 3 a"} which is the desired NTB. The synthesis, however, avoids grammar con¬ 
struction as well as specification of the language that guarantees the desired regenerative 
behaviour, e.g. pr{B), which is necessary for other approaches. The parameter k here, in 
fact, parametrises the family of non-blocking regenerative supervisors for the system in 
the example given. 
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Abstract. This paper is in continuation of our earlier work on the role of 
hydrodynamic stability theory in understanding wall-bounded turbulent flows. 
Work in this area was pioneered by Malkus, followed by Reynolds, Tiederman 
and Hussain. Numerical results showed that the linear instability inodes are 
damped, a result also confirmed by our earlier work for the boundary layer 
flow case. This led to waning of interest in this approach. 

In the present work the problem is reformulated using an improved non¬ 
isotropic model for the stress tensor based on the model of Pope. This improved 
model does yield unstable wall modes. A wide range of unstable wavenum¬ 
bers is observed and these unstable modes mimic some of the key features of 
wall-bounded turbulent flows. Comparisons with experimental data are also 
presented. The present work keeps alive the question of relevance of stability 
theory to wall-bounded turbulent flows. 

Keywords. Hydrodynamic stability; wall-bounded turbulent flow; boundary 
layer flow; non-isotropic model. 


. Introduction 

'he motivation for the present paper is primarily to re-examine a question that has been 
sked before, but has got mixed answers so far. Whether or not hydrodynamic stability 
leory has any bearing on the problem of wall-bounded turbulent shear flows, especially 
2 view of the fact that it does have a significant role in free turbulent shear flows. See, 
or example, the discussion by Caster et al (1985) on the turbulent mixing layer, the 
omprehensive review on stability and free turbulent shear flows by Liu (1988), and a 
lore recent overview by Roshko (1992). 

The question of a connection between stability theory and turbulent shear flow was 
irst raised by Landau (1944) based on a nonlinear stability model. While his work did 
lOt prove to be a suitable model for turbulence, his equation for nonlinear growth found 
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many applications in the field of instability and transition. The next work of great concep¬ 
tual importance is that of Malkus (1956), wherein a theory of turbulence was developed 
based on the introduced concept of marginal stability. According to this concept, it was 
proposed that if the mean velocity profile typical of wall-bounded turbulent flows is 
used in the solution of the classical Orr-Sommerfeld equation, then the profile would 
prove to be marginally or neutrally stable at the existing flow Reynolds number. It is 
significant also that Malkus suggested that the molecular viscosity and not the eddy vis¬ 
cosity be used in the solution of the Orr-Sommerfeld equation. Malkus’s theory was 
rigorously put to test in an important work by Reynolds & Tiederman (1967). This work 
also gives a lucid review of Malkus’s theory. Reynolds & Tiederman (1967) investi¬ 
gated the stability of fully developed turbulent flow between parallel plates, on the lines 
of Malkus’s proposed theory. They used the turbulent mean velocity profile for chan¬ 
nel flow in the solution of the classical Orr-Sommerfeld equation (using the molecular 
viscosity in the Orr-Sommerfeld equation). The results obtained showed without doubt 
that Malkus’s theory as proposed was not valid, and there was a huge discrepancy be¬ 
tween the flow Reynolds number and the Reynolds number corresponding to neutral or 
marginal stability, the latter being an order of magnitude higher than the former. Prabhu 
(1968, private communication) also independently obtained the same result. More re¬ 
cently, Sen et al (1993) obtained the same result for the turbulent boundary layer flow 
problem. 

Although Malkus’s theory was not proven, it did not immediately detract attention from 
the general question of connection between instability and turbulence, and, as mentioned 
before, considerable success was obtained later on in linking the two for the cases of free 
turbulent shear flows. The results for free turbulent shear flows were very decisive mainly 
because of the inflectional mean-velocity profile present in such flows. Thus the dominant 
instability is the inviscid instability, which continues to exist even in the presence of the 
existing turbulent flow field. Also, the dominant eddy structure of the turbulent flow field is 
found to correlate well with the vorticular structure corresponding to the inviscid instability 
(see Liu 1988 and Roshko 1992 for comprehensive reviews). 

The general question of connection between instability and turbulence, with reference 
to wall-bounded turbulent shear flow, went through another round of serious examination 
by Reynolds & Hussain (1972). This time, abandoning Malkusian precepts, first the basic 
equations were derived for a superposed organised disturbance in turbulent flow. Thereafter 
an equation for the stability of this organised disturbance was derived, which was like 
an augmented form of the classical Orr-Sommerfeld equation, but containing extension 
terms dependent on the eddy viscosity. Underlying the model was a closure problem for the 
Reynolds stress tensor, which was resolved based on what the authors called the ‘Newtonian 
eddy’ model. Results of solution of their extended Orr-Sommerfeld equation, for channel 
flow, again yielded damped modes. Nevertheless, the results were closer to establishing a 
connection between instability and turbulent wall flows than was obtained in Reynolds & 
Tiederman (1967). Experiments performed by Hussain & Reynolds (1970,1972) showed 
some agreement with the theory of Reynolds & Hussain (1972). Subsequent to the works 
of Reynolds & Hussain (1972), some nonlinear and three-dimensional theories were also 
developed (see for example the work of Jang et al 1986), but these are outside the scope 
of the present discussion. 



The present work seeks to re-examine the question of connection between instability 
nd wall turbulence. The organised disturbances are linear and two-dimensional, and the 
articular flow chosen is the turbulent boundary layer flow over a flat plate. The chief 
nprovement over the work of Reynolds & Hussain (1972) is that a more realistic and 
nproved model is chosen for the turbulent stress tensor, based on the anisotropic model 
f Pope (1975). This model gives further extensions of the Orr-Sommerfeld equation 
ver what Reynolds & Hussain (1972) obtained. The results show that an unstable wall- 
lode exists over a wide range of the spatial wavenumber a. The instability characteristics 
3ale very well with the inner variables of turbulent flow. Also the organised disturbances 
limic some of the key features of wall-turbulence. Experiments performed for the one 
imensional energy-spectrum for the turbulent longitudinal velocity, indicate that the range 
f unstable wavelengths is well contained within the energy-containing part of the energy 
Dectrum. Thus, the present work keeps alive the question of relevance of stability theory 
1 understanding turbulence in wall-bounded turbulent flows. 

Preliminary versions of the present work appeared in Sen & Veeravalli (1994, 1996). 


. Theory 


’he instantaneous velocity vector w, obeys the Navier-Stokes and continuity equations: 
Sui dui 1 dp d^Ui 

— u j — =-4- y--—, (1 < 

dt dxj p dxi ^Xjdxj 


dui 

dXi 


= 0 . 


(lb) 


’he velocity and pressure fields are decomposed in turbulent flows by the well known 
Leynolds decomposition; typically, 

Ui=Ui+u[-, p=T + p'- (2) 

lere Ui,~p are respectively the mean velocity and pressure, and u\, p' are the (random) 
Luctuations. If we now superpose an organised (solenoidal) disturbance ui, p (with zero 
nean), the instantaneous velocity and pressure are respectively given as follows: 

ui =Ui-\rUip=~p + p + p. (3a) 


"he time averages of ui, p are still respectively m,-, p, but the ensemble (phase-locked) 
tverages are different. The ensemble averages are given respectively as 

{ui)=Ui+Ui-, (p) = p + p- (3b) 

d[oreover, the organised disturbance is assumed to be small, or linear, and in addition it 
)beys the following: 

\{uiUj)\ <\{u'iu'j)\. (4) 

rhe above assumption (4) restricts the organised disturbances to being weaker than was 
considered by Reynolds & Hussain (1972). 
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For future clarity some definitions and notations are introduced: (i) an overbar (-) over 
any quantity implies time average; (ii) the symbols () enclosing a quantity denote ensemble 
average. Further, we note that: (i) when the flow field is not perturbed by an organised 
disturbance, the time average and the ensemble average are the same; (ii) when the flow 
field contains an organised disturbance, the time average and the ensemble average are 
different; (iii) it will be seen later that if the organised disturbance in the flow field is small 
(in the sense defined in (4) above), then the time average is the same as the time average of 
the con'esponding flow field that is not perturbed by the organised disturbance. Hereafter, 
the terms ‘perturbed flow’ and ‘unperturbed flow’ will mean respectively the flow that is 
perturbed and that is not perturbed, by the organised disturbance. 

One needs now to look into the forms that (la) takes upon different ways of averaging. 
First the typical time average, with the disturbance decomposition as in (2), which yields 
the well-known Reynolds averaged equation: 

dui _ dui 1 

--f- u j -—-“t” n-—-. (5) 

dt dxj p dxi dxjdxj dxj 

One can also obtain the ensemble average of (la), to obtain the equation for (w/). We need 
to remember the following rules for ensemble averaging, given with respect to any two 
generic quantities a, b. 

(db) =db; {a'b) = 0; (6) 

and, in view of (4), and in the event that d,b 0{ui), then the following is also valid: 

{db) 0; db ^ 0. (7) 


In view of (4) and (7), it may be noted that the time average of (la), even when distur¬ 
bance form (3a) is being considered, becomes the same as (5). Particularly, the organised 
disturbances are assumed small enough not to alter u'-u'^. Further discussion on this point 
will be taken up later. 

Next we look at the form of (la) after ensemble averaging, with the disturbance decom¬ 
position as in (3a). Remembering (6) and (7) one obtains: 


9(«/) , , ,9(wi: 

—— -h(Mj)-— 

at dxi 


1 dip) d^iui) 
-—— + v 
p axi 


diiu'.u'.)) 


dxj dxj 


dXi 


( 8 ) 


One is now in a position to obtain the dynamic equation for the organised disturbance, by 
subtracting (5) from (8), and remembering (3b). This gives the following equation: 

B^ui dra 

■ ' (9) 


dui __ dui 
dt dXj 


„ dui \ dp 

“f Uj - - =--b V 


dx 


.j p dxi 


dxjdxj 


+ 


dxj ’ 


where r/j is the difference between the ensemble average and time average of u^-Uj, that is 
it is the modulation in the Reynolds stress tensor, and is given as follows: 




j. It will be instructive to look at simplified physical models and reasoning to understand 
e manner of resolution proposed for the closure problems. 

At this stage we introduce (twice) the rate of strain tensor and (twice) the vorticity tensor, 
spectively stj and ojij , (in generic form) as follows; 


s 




( dui 

dxj dxi j ’ 


0)1 r = 


( dui 



(11 a,b) 


he expressions in (11) are in generic form, implying that if, for example, u is replaced 
/ u, then Si j and coij are respectively replaced by Sij and Wij. 

Let us now restrict our attention to 2-D near-parallel flows, specifically the turbulent 
3undary layer. A definition sketch is shown in figure 1. Note that in the discussions to 
>llow, the vectors (u\,U 2 ,ii 3 ) and {u,v,w), will be used interchangably and so also 
:\,X 2 ,X 2 ,) and (x, y, z). The x coordinate is in the direction of the free stream, y is 
L the direction normal to the wall, and z is the transverse direction. Also we make the 
jasi-parallel assumption due to which the mean-velocity field is given as if = w(y), 
0, and uJ = 0. Further, all the mean-velocity gradients, except for du/dy, are either 
iro, or negligibly small. The turbulence (M-)-field is also assumed homogeneous in the 
direction, and near-homogeneous in the x direction. All derivatives of time-averaged 
Liantities are zero in the z-direction, and nearly zero in the x-direction. Moreover the 
relations u'w^ and v'w' are zero. However, is non-zero, and this term keeps the 
eynolds stress tensor u'^u'- three dimensional. The velocity scale for the problem could 
e either the free stream velocity, Voo (outer scaling), or the friction velocity, u* (inner 
mling). The relevant length scales are the boundary layer thickness 8 or the displacement 
dckness 8* (both outer scaling), or the viscous length scale, v/u,,, (inner scaling). The 
laracteristic Reynolds number for the problem is R = Uoo8/v. 

With the problem posed as above, two closure schemes are presented. The first, called 
le ‘isotropic eddy viscosity model’ is identical in form to the Newtonian eddy model of 
eynolds & Hussain (1972). The second, called the ‘anisotropic eddy viscosity model’ is 
ased on the model proposed by Pope (1975). 


2.1 The isotropic eddy viscosity model 


The isotropic eddy viscosity model, for the unperturbed flow, is given by the following 
well-known expression: 

—u'iu'j = —\k8ij + with, k = (12a,b) 

In (12), k is the turbulent kinetic energy, and 6 is the eddy viscosity. Also, in the wall- 
bounded two-dimensional flow being considered here, € = e{y), that is, the eddy viscosity 
is a function of the coordinate perpendicular to the wall. Details of the form of € will be 
discussed later. 

Equation (12), despite similarity in form with the corresponding form for laminar sheai’ 
flow, is no more believed to represent a gradient-type diffusion of momentum. In fact, the 
justification for the above form is the following (see Tennekes & Lumley 1972).-First, the 
eddy viscosity is of the form e ~ u'’/, where I is a turbulent length scale (which may loosely 
be called the mixing length for historical reasons). This length scale is unique for wall- 
bounded turbulent shear flows, like in the case of the flat plate, and is given as / ~ (/cy) in 
the inner region and approaching a constant in the outer region. Also k is the Von Karman 
constant. Second, the turbulent velocity scale is given uniquely by u* = where 

Tyj is the wall shear stress; hence u' scales with Third, the turbulent time scale and the 
mean flow time scale are of the same order. Thus, form (12) emerges as a dimensional 
necessity. Since it is not the scope of the present paper to go into a detailed discussion on 
these points, it suffices to say that form (12) is a valid and legitimate form for wall-bounded 
shear flows, like the flat plate, that are governed by one velocity scale and one length 
scale /. 

The ideas developed over the past years, on lines of the above discussion, can be sum¬ 
marised by the proposition of the ‘generalised eddy viscosity hypothesis’ (GEVH). Such 
cases of turbulent shearflow that have one characteristic velocity scale and one character¬ 
istic length scale of turbulence and, where the turbulence is sustained by the mean shear, 
will have a time scale of turbulence that is of the same order as the mean flow time scale, 
the latter being given by the reciprocal of the dominant mean-velocity gradient. Further, 
the turbulent shear stress tensor in such flows is capable of being described by products 
of scalar functions derived from the turbulence and appropriate tensorial combinations 
of the mean-velocity gradients. Further, any description of the shear stress tensor must be 
tensorially consistent, that is, it must obey the laws of tensor transformation. In actuality, 
form (12) is the simplest and most widely used manifestation of GEVH. 

What is important for us to note, in the context of the present work, is that if (12) be 
looked upon as a constitutive equation for u\uj, then the right hand side of this equation 
has one part that is solely dependent on the existing turbulence, viz. € ~ v'l, and the 
other part on gradients of the mean flow. It is reasonable to assume that if such a flow is 
perturbed, yet leaving u' and I unchanged, then the eddy viscosity c remains unchanged. 
In other words, if the existing turbulence is unchanged then c remains unchanged. Further, 
if the perturbations cause modulations in the mean field, then these modulations appear 
as modulations of the Reynolds stress tensor. The organised disturbance does exactly this. 
First it modulates the mean flow, as may be seen from (3b). Second, it can be seen from 
Reynolds & Hussain (1972) that the modification of the turbulent kinetic energy, due to the 



presence of the organised disturbance, is 0{u^). In view of (4) this is negligible, and thus 
there is no change in the turbulence, and consequently no change in the eddy viscosity. 

In the extreme case of a very high-amplitude organised disturbance, possibly even the 
form (12) would not be valid, since the imposed disturbance would introduce its own length 
and velocity scales in the problem, thus negating the validity of GEVH. 

The consequence of the above discussion is significant, in that (12a) can be extended to 
the flow modulated by the organised perturbation in the following manner: 

- {u'iu'j ) = ^ e{sij ). ( 13 ) 

It is easy to see now that (u'-Uj) will be modulated by the organised perturbation despite 
e having remained unchanged. It is assumed in (13), following Reynolds & Hussain (1972) 
that the organised disturbances do not modulate the dilatational part of the stress tensor. 
Thus {{k) — k) is neglected, where {k) = (1 12){u'ju'j). Besides, even if ({k) — k) were to 
be retained, it would be like the the pressure term in (9), and this would not eventually 
affect the derivation of the extended Orr-Sommerfeld equation. Further, by subtracting 
(12a) from (13), it is easy to see that rij ~ 0{ui), which is consistent with Reynolds & 
Hussain (1972), and the expression for rij is given as follows: 

fij = €Sij. (14) 

Reynolds & Hussain (1972) have not considered m, as small as in (4) herein, nor have 
they considered very high amplitude perturbations. However, they have reasoned that, 
with some restrictions, in their problem also the Newtonian eddy model will give the same 
expression for fij as in (14). 

2.2 The anisotropic eddy viscosity model 

The necessity of considering an improved model, over the isotropic eddy viscosity model, 
is quite strong, and various angles of the matter will be looked into. The main reason for 
proposing an improved model (without keeping any allied problem, like the hydrodynamic 
stability problem, in view) has been given by Pope (1975). There could be some flows, 
wherein the turbulence is isotropic or anisotropic, where GEVH is not valid (one such 
simple example is given by Tennekes & Lumley (1972), of turbulent flow over a flat-plate 
with vertical blowing or suction through the plate). Also, there are some flows wherein 
the turbulence is anisotropic, yet GEVH is valid. The most obvious examples in the latter 
category are the simple wall turbulent flows, the flat plate boundary layer problem being 
one of them. The important contribution of Pope (1975) has been to develop advanced 
level constitutive equations under GEVH, mainly to account for anisotropy in flows that 
obey GEVH. 

First we examine the shortcomings of (12) as a constitutive equation for the Reynolds 
stress tensor. It can be seen from (12a) that the normal stresses predicted by this model 
give equal to ^k. This is quite a departure from the known anisotropic 

behaviour in the entire shear layer and particularly in the wall region. To overcome this 
shortcoming, the constitutive equation that has been given by Pope (1975), under the aegis 
of GEVH, is as follows: 

-u'iu'j = -\k8ij -h €Sij - €C{kj€d)[\\a)ikSkj 


(15) 
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In (15), €ci is the dissipation rate of the turbulent kinetic energy, and C is a constant to 
match experimental results. The last term in (15) accentuates the anisotropy generally by 
augmenting reducing and leaving unchanged. Further, the last term leaves 
u'v’ unchanged. Pope (1975) has one more term in (15), involving a scalar invariant of 
Sij, namely the scalar product 'sijSji. This term tends to equalise and v''^, but makes 
w''^ different from both. This term is rejected for the disturbance equation on two grounds. 
First, the accentuation of anisotropy in the z direction has not much relevance in the 
present problem in which basically two-dimensional disturbances are being considered. 
The second reason is more serious. The SijSji-term has to be incorporated artificially with a 
two-dimensional Ki'onecker delta. Since the latter does not behave like the identity matrix 
in three dimensions, it destroys the tensorial symmetry of the problem, and, eventually 
renders the disturbances non-transformable by Squire’s theorem. Pope (1975) himself has 
suggested that this term can be ignored if accentuation of anisotropy in the z direction is 
not important. Thus, in tlie present problem, the SijSji-t&rm is nominally considered in 
the (unperturbed) mean-motion only, mainly to estimate the value of the two anisotropy 
parameters C and Ci, the latter to be introduced later. 

Equation (15) has a unique tensorial form involving T/y and Wij under GEVH. Over 
and above the basic isotropic eddy viscosity model, (12), this is the next advanced-level 
constitutive equation, from amongst the hierarchy of constitutive equations derived by 
Pope (1975). So long as GEVH is valid, and a model is to be proposed within GEVH, the 
uniqueness of the form of (15) leaves little scope for seeking out other tensorial forms at 
the same advanced-level as in (15). At most, one can look for alternative but equivalent 
forms for the turbulence-derived scalar functions like k/cd- It is therefore reasonable to 
conclude, that, since the present problem obeys GEVH, form (15) is clearly both a unique 
and logical improvement over (12). Also, other such turbulence models that are not within 
the framework of GEVH, are outside the scope of the present paper. But then, such models 
would not also admit even a form like (12), and it is doubtful whether such models would 
capture the real physics in the present problem, as much as GEVH would. 

Whilst retaining the basic tensorial form as in (15), for certain conveniences in calcula¬ 
tions we introduce an alternative but equivalent scalar function in place of the turbulence- 
derived scalar function k/^d in (15). Now, k/cd represents a turbulent time scale. It is 
convenient to put this term equivalently as k/cd = A/(m)', where («)' = (du/dy). The 
function X is called the anisotropy function, with the constant C in (15) absorbed in A.. 
The anisotropy function A thus physically represents the ratio of the turbulent time scale 
to the mean flow time scale. The constitutive equation adopted in the present problem is 
therefore given as follows: 

-u'iu'j = -^kSij -1- esij - €(X/u')[^[djikSfcj - Sikcdicj]]. (16) 

Using exactly the same arguments as in § 2.1 earlier, (16) may be extended to perturbed 
flow, as follows: 

-(u'^Uj) = -\k^ij + - €(X/u'){\[{a)ik){skj) - {sik){oikj)^]■ (17) 

Also, as before, rij may be obtained by subtracting (16) from (17), and is given as follows: 



2.3 The anisotropy function 


A word is now in order regarding the anisotropy function A, and its efficacy in the present 
formulation. One can put ed ~ u’v'(u)' in (16), the expression for €d becoming almost 
exact in the outer region where production nearly equals dissipitation. The expression for 
A is thereafter given as follows; 

A = C(k/(~^)), (19) 

where C is either a constant, or at most a slowly varying function. It is noted from (19) that 
A is a universal non-dimensional function, and is like the inverse of the structure function 
given in Townsend (1976), and also plotted by Hinze (1975; figures 7.21 & 7.22). In fact, 
sin ce k /ed in (15) is meant to extract a turbulent time scale, one might as well replace Cd 
by u'v'(U)' in the entire range of y in (15). 

The contribution of the 5/yTy/-term, to the (unperturbed) mean-motion, is nominally 
added to (16) using a form similar to (19) with C\ in place of C. Thereafter, upon expanding 
(16), one obtains the following expressions for the turbulence stresses: 
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The accentuation of anisotropy now becomes clear, and 
be rewritten as follows: 


i/2 


,/2 


= - + C-4Ci; — = --C4-Ci; 
k 3 k 3 


(20a) 

(20b) 

(20c) 

(20d) 

in view of (20d), (20a,b,c) may 
w'^ 2 

= --2Ci. (21a,b,c) 

k 3 


Estimation of the values of C and C\ have been carried out from Klebanoff’s experimental 
data (cf. Hinze 1975), and from the DNS data of Spalart (1988). Thereafter, using (19), A 
has been evaluated. The results are shown in figure 2a, in which y+ = yu,^/v is y-scaled in 
inner variables. It is seen that Ci is quite small except very close to the wall, where it has a 
value of ~ 0.2. In case that Ci = 0, then it is seen from (21a, b, c) that there is a restriction 
on admissible values of C, viz. C < 0.67, to avoid negative values of v'^. However, close 
to the wall, where Ci ~ 0.2, C can have values ~ 0.9, as may be seen in figure 2a. 

It is seen in figure 2a that A has a value of 25 close to the wall, relaxing to around 2 in 
the region y"*" ~ 60. In the outer region it relaxes further to around 1, and drops to 0 close 
to the edge of the boundary layer where there is isotropy in the sense that w'^, and w''^ 
are nearly equal to each other. The relatively high values of A obtained close to the wall 
are because of the strong anisotropy in the wall region. 

Stability calculations were performed with a wide latitude of variation in the form of 
A, retaining the feature that A has a high value in the wall region. It was observed that 





(b) y 


Figure 2. (a) Graph of anisotropy function A"''(= A.) versus standardised in terms of 

inner variables. Based on experimental data of Klebanoff, and from (22). Also shown are C 
and Cl. (b) Graph of anisotropy function X versus y{= y&/8), standardised in terms of outer 
variables, based on (23). 
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the exact shape of the A,-curve is not important. Rather, instability is obtained with a high 
value of X (~ 9.0 or more) in the wall region. Thus instability is due to the anisotropy in 
the wall region, a feature that could not have been captured by the isotropic eddy viscosity 
model. 

Two of the forms used for X in the calculations are given below. The first is based on 
inner variables scaling and is given as follows: 


= >.h 




1 Xlo/Xfii 


1 — erf 


Ox 


with. 


"Xhi — 10.5; Xio — 3.0; — 18.0, o\ 7.5 


( 22 a) 

( 22 b) 


The A.-curve based on (22) is shown in figure 2a. Equation (22) renders the stability equa¬ 
tions stiff and the calculations are more laborious and constrained. This is because the 
region of high X is squashed close to the v/all. For facility in calculation another expres¬ 
sion for X was used based on outer variable scaling (plot shown in figure 2b) which gives a 
value of 10 at the wall relaxing much more gently to a value of 3 in the outer region. This 
equation is given as follows: 


^(3^) = A./I/ 
with, 


^loIXhi + 


1 ^lo/^hi 


)l'- 
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Xhi = 10.0; Xio = 3.0; yox = 0.5; = 0.25. 


(23a) 

(23b) 


Otlier forms were also tried, varying the parameters in (23). The fact that instability is not 
dependent on one particular critical choice of the shape of the A,-curve, but is generally 
governed by high values of X in the wall region, is a point that greatly favours the statement 
of a general proposition: that instability is due to the anisotropy in the wall region. 


2.4 The disturbance equations 

We simplify the evolution equation (9) for «/, remembering that the disturbances are two- 
dimensional and Squire-transfonnable, and that the mean-flow is amenable to the quasi¬ 
parallel approximation. A stream function f is introduced for the organised disturbances 
such that u = dxlz/dy and v = -9 i/f/3x. After assuming normal modes, i/r may be 
expressed in the following form: 

iff = 0 ( 3 ;) &xp[ia{x — ct)], (24) 

where, a is the spatial wavenumber and c = Cr-{-ici is the (complex) wave speed. Introduc¬ 
ing (14) in the evolution equation ( 9 ) for ui , remembering the quasi-parallel approximation 
for the mean flow, and remembering the closure equations (14, 18) for rtj, one arrives at 
extended forms of the Orr-Sommerfeld equation (a more detailed derivation may be found 
in Sen & Veeravalli 1994). For the isotropic eddy viscosity model this is given as follows: 

ial(u - c)(4>" - oP-<t>) - «"<#>] - ( 1 /R)[ 0 "" - 2 oV" + «'*'!(>] 

-(\IR)[E[<j>"" - 2a^4>" + + 2£'|0"' - cp-<l>'\ 


follows: 


ia[iu - c)((l)" - a^(l)) - w'Vl - - 2 q;^ 0 " + a^</)] 

-\/R[E{(l)"" - lop-ip" + + 2E'{(j)"' - a'^cf)'} 

+E"{(j)" + a'^cj)}] + (XE/R)[-2ici(l)"' + 2ia^(t)'] 

. -j-(2ia(p'/R)[XE'' + 2X'E' + X"E] = 0. (26) 

In both (25) and (26), primes (') denote differentiation with respect to 3 ;. All quantities 
in (25) and (26) have been non-dimensionalised by C/ 00 , <5 and y, i.e. E = €/v,u = 
Wd/ C/ 00 , y = yd/5 etc., where subscript ‘d’ represents dimensional or physical quantities. 
Also, in (25) and (26), the first group of terms in square brackets corresponds to the 
Rayleigh equation; the first two groups of terms in square brackets correspond to the 
classical Orr-Sommerfeld equation; and the remaining terms constitute the modifications 
in the classical Orr-Sommerfeld equation. Now, E” -> 0 (to be shown later) at the edge 
of the boundary layer, and, X can be made to obey the less stringent conditions X', X" —> 0 
(without letting A, -» 0) at the edge of the boundary layer. Thus, the appropriate boundary 
conditions for the two extended Orr-Sommerfeld equations, viz. (25) and (26), are the 
same as those for the classical Orr-Sommerfeld equation, namely: 

0 ' = 0 at y = 0; (27a) 

0(y) ~ for y -> CO. (27b) 

Either of the disturbance equations (25) or (26), along with the appropriate boundary 
conditions (27), constitutes an eigenvalue problem. In the temporal problem, a and R 
are chosen real and c = c,- + ici is determined as the (complex) eigenvalue. Stability or 
instability is obtained respectively as Cj < 0 or Ci > 0. As noted earlier the extended 
Orr-Sommerfeld equation, as given by (25), did not yield any unstable modes although 
near neutral modes were observed. 

Next we look at the inner-variable form of (26). It is well known that at high values of 
R, the turbulence in the region close to the wall assumes a universal character, which is 
best seen by expressing the governing equation in inner variables. The velocity and length 
scales are and y/unc, respectively. Denoting quantities scaled by inner variables by the 
superscript (+), we note that, X = A+; E — B = vl/U^; y"^ = yR-/B-,u'^ = 
u/v^; and, = a/{R-/B), where B is the non-dimensional wall shear stress. Thus, in 
inner variables, the disturbance equation (26) becomes: 

/a+[(M+ - c+)(0" - Q!^-^0) - If+'V] - W" - 2 q!+V" -b 

-[E"^{<1)’"' - 2 q;+V" + 2 £'+'{ 0 '" - Qf+^ 0 '} 

4 . + A+E+[- 2 /q'+(/)'" + 2 ia+V'] 

-|-2/q!+</>'[A.+E+" + 2A+'e+' -f A+"e+] = 0; (28) 

where, in (28), primes (') denote differentiation with respect to y"^. Note that (28) does 
not show any explicit dependence on the Reynolds number, R. As will be seen later, (28) 



corresponding equation (28) would indicate that the behaviour of these wall modes should 
be near universal at high R. This is because the outer region does not have a strong influence 
on wall modes. Universality is actually confirmed by the present numerical results. 


3. Details of various functions 

The mean flow in wall turbulence is known to be multi-layered. Closest to the wall is 
the viscous sub-layer, followed by the buffer zone, and then followed respectively by 
the inner and outer log-law regions. Mathematical expressions are available for the mean 
velocity distributions, individually for each of these respective zones; but taken together 
these expressions do not constitute an analytically continuous expression for the mean 
velocity distribution. For hydrodynamic stability calculations as in the present problem, it 
is necessary to obtain an analytically continuous expression for the mean flow. 

One very successful method of generating an analytically continuous mean velocity 
distribution for turbulent channel flow is reported in the work of Reynolds & Tiederman 
(1967). In the present work, the mean velocity profile for turbulent boundary layer flow 
is generated on lines of the work of Reynolds & Tiederman (1967). This requires obtain¬ 
ing suitable expressions for the eddy viscosity, and the shear stress distribution. This is 
discussed next. 


3.1 The expression for eddy-viscosity 


First we look at the eddy viscosity E(y). The expression for the eddy viscosity E{y) is 
given by a combined expression of Van Driest’s law for the wall region and Reichardt’s 
expression for the outer region. The composite E{y) is given in Reynolds & Tiederman 
(1967), and Reynolds & Hussain (1972). Also, a slightly different expression is given in 
Hussain & Reynolds (1975). The Reynolds & Tiederman (1967) expression is given as 
follows: 




1 + K^^-^(2y - y'^fo ~4y-\- 2y^)^ 


X 11 — exp 


-yRVB \ 
A+ j 


(29) 


In (29), K is the Von Karman constant, A+ is the Van Driest constant. The above equation 
is best suited for channel flow and a detailed plot of E (y) may be seen in Reynolds & 
Tiederman (1967). A slightly modified basic expression for £(y) is needed for boundary 
layer flow, because E', E" ~ 0 for y = 1 (boundary layer edge). This basic expression, 
denoted by Eb, again needs to be further modified to account for the intermittency in the 
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outer region. But first we define the basic expre ion Eb, as follows: 


Eb(y)=:^ 


1 + at" 
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■(2y-y^) (2-2y + y^y 
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1 1 - exp 


(^) 
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2 ‘ 


( 30 ) 


In the expression for i.e. in (30), k has been chosen as /r = 0.4, and A"*' has been 
chosen as A^" = 27. Next we see in (30), what further modifications in Eb are needed 
in order to accommodate the intermittency in the outer region. Although Eb satisfies 
E'jj, E'^ ~ 0, at y = 1, it still does not approach the right numerical value of £■ at y = 1. 
So, the first correction is to multiply Eb by an error function g(y), so that the right value 
of the eddy viscosity is obtained at y = 1, prior to correction for intermittency. Next, to 
account for external intermittency of the boundary layer, for y > 0.4, the eddy viscosity 
expression is further multiplied by the intermittency function (another error function) y (y). 
The final expression for eddy viscosity is given as follows: 

E(y) - Eb(y)8(y)yiy)', (31a) 


where, 
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with, a = 0.2; yog = 0.4; ag = 0.13; 
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and, 
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yoy = 0.797; 

Plots of E (y), for different values of R, and corresponding to (31), are shown in figure 3a. 
Also, figure 3b shows the plot for e/(u*(5) for R = 75, 000, and, shown in this figure also 
are the points coiTesponding to the measurements by Klebanoff (1954), reported in Hinze 
(1975). 

One may look at the form of Eb in inner variables, to make an assessment of the 
(marginal) departure from universality, of the disturbance equation (26) and (28) because 
of the outer scaling. Thus using inner scaling in (30) we have 
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Note in (32) that y = y+/(/?\/5). Therefore, the outer region does contribute to E^, 
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Figure 3. (a) Graph of eddy viscosity E(= e/v) versus y(= yi/S). (b) Graph of eddy 

viscosity €/(u*5) versus y(= yd/<5). Comparison with Klebanoff’s data at Re = 75,000. (c) 

Graph of eddy viscosity £■*■(= E) versus y"*", inner variables. Scale of y+ is logarithmic. 

wall mode that is obtained in the problem. One may now similarly express (31) in inner 
variables scaling, to obtain (which is the same as £’(y’^)). Figure 3c shows plots 

of the eddy viscosity E'^ (y"*"), in inner variables. This figure shows that in the inner region 
the eddy viscosity is the same across different R, a feature that is responsible for the 
universality of the wall modes. 

3.2 The expression for the shear stress distribution 


The distribution for the (total) shear stress r is needed in order to be able to integrate the 
equation of motion to obtain the mean velocity distribution u(y). It is well known from 
experimental data that the distribution of mean shear looks similar to a tanh-profile, or a 
complementary error function profile. The expression chosen here is again a composite 
expression comprising an error function and a Gaussian function, to be able to match the ex¬ 
perimental shear stress distribution, and, also, for the mean velocity profile to be able to sat¬ 
isfy the proper outer boundary conditions. This composite expression is given as follows: 

r(y) = rujFiy)Giy); (33a) 


where, Tm is the wall shear stress. Also r is dimensional or non-dimensional depending on 
whether z-uj is dimensional or non-dimensional. When Xyj is non-dimensional then Zy^ = B. 
Further, 


F(y) = l+b 


1 -f- erf 


y ~ yoF j 


gf 


with. 


b = 0.7; yo/r = 0.56; ap = 0.204; 


(33b) 






Figure 4. Graph of shear-stress t/tu, versus y„i. Here y is rescaled so that at y,,, = 1, 
u = 0.99. Data points of Klebanoff are also shown, for Re = 75,000. 


and, 


G();) = 


1,G'(3;)=0; forO 


exp 




3^0G J 


< yoG‘, 

for y >yoG\ 


with, Jog = 0.05; og — 0.26. 


(33c) 


A typical plot of v/tu] versus j is shown in figure 4. In this figure the scale for j is modified 
to jm (defined in the next sub-section). Also shown are data points for r/xyu measured by 
Klebanoff (1954), and as reported in Hinze (1975). 


3.3 Generation of the mean velocity profile 


The x-momentum equation for mean velocity simplifies for boundary layer flow to the 
following: 


_du _ 

u -h V 

dx 


du 

9y 


iA 

^aj 


a + E) 


du' 

9y 


dr 


(34) 


In the discussion that follows, the convective acceleration terms are not required. Hence 
(34) gives a relationship between E and r (which is in non-dimensional form). Thus, after 
integrating (at any given station in x), (34) becomes the following: 


du 

a + E)—=:Rr-, 
dy 


(35) 




Figure 5. Graph of mean-velocity u versus . y rescaled so that at y,„ 
Data points are from present experiments, for Re around 17,000. 










where, in (35), r is given from (33a,b,c). Also, it may be noted in (35) that du/dy satisfies 
the condition du/dy = RB,for y = 0. Also du/dy satisfies the condition du/dy 0 for 
y -> oo; because, £ 0 for y -> oo, and, the right hand side of (35), viz. Rr 0 for 
y oo. Integrating (35) we obtain 


RB 


L 


y F(y')G(y) 
0 1 + E(y') 


dy'; 


(36) 


where, in (36), F(y), G(y) are given from (33a,b,c). Since u must approach the free stream 
velocity Uqq (= 1 in non-dimensional form), therefore, it follows from (36) that, 


1 = RB 



F(y)G(y) 
1 + E(y) 


(37) 


Thus iterating between (30), (3 la,b,c), (33a,b,c), and (37) one obtains the converged value 
of B. After this the expression for u is given from equation (36). 

Typical plots of the mean velocity distribution u, both in outer and inner variables 
are shown respectively in figures 5 and 6. Also shown in these figures are the points 
corresponding to the present experimental measurements. Solution for u from (36) is 
given in terms of y(= yd/^)- Thus y = 1 represents yd = 5. Now, after w has been 
calculated from (36), it is not guaranteed that at y = 1, « has the exact value oiu = 0.99. 
Therefore, for comparison with experiments, it is sometimes expedient to rescale y to y,„, 
where y^ = y/^^ 99 , and, for y = d?99, Ti = 0.99. The scale for y, in figures 5 and 6, is 
thus modified to y,n, for ease of comparison with experimental results. Since the Reynolds 
numbers (in figures 5 and 6) are relatively small, the inner and outer regions are not well 
separated. Therefore, one need not expect the log-law to be valid over an extended region 
in y^. This may also be confirmed from the DNS data of Spalart (see figure 5 of Spalart 
1988). In figure 6 it may be noted that there is very good agreement in the inner-wall 
region. This is what is really important in the stability calculations because the critical 
point lies in the inner-wall region. In the outer region also, there is good agreement with 
the experimental data points. 

It is clearly seen from the above figures that a very satisfactory method of generating an 
analytically continuous mean velocity has been obtained. Further refinements, if needed, 
may be made by adjusting the constants in the expressions for E and t. 


4. Numerical results 

The numerical solutions of the extended Orr-Sommerfeld equation (26) were based on 
finite-difference methods developed in the earlier works of one of the authors (see Sen 
& Vashist (1989) and Sen 1993). In the present problem, the methods described in the 
earlier works for the stability of the Blasius boundary layer are extended to the present 
turbulent boundary layer problem. The resolution used is ten times greater, and the basic 
step size in y, in the present problem, is /i = 0.001. Over a range of 0 < y < 1.5, 1501 
points are required for each array specification. A seven-point finite difference scheme, 
employing a Noumerov transform and using a molecule described in Sen et al (1985), 
was employed in the calculations, and this gives basic round-off accuracy as 0{h^), i.e. 



Figure 7. Graph of eigenfunction 0 = + i0,- versus y. Normalisation (j) — 1 + Of at 

y = 1.5. Also Re = 5000, ol = 2.7, Cr = 0.3486034, q = 6.0044436. 









Figure 9. Graph of r/? du/6y versus y. Normalisation, (j) = 1 + Oi at _y = 1.5. Also, 
Re = 5000, a = 2.7, c, = 0.3486034, q = 0.0044436. 


10“ with h = 0.001, Thus all care was taken to ensure that the very rapid changes taking 
place in the inner region, are faithfully taken into account in the numerical calculations. 
Also, the Blasius boundary layer results were verified as a cross-check. Calculations were 
performed using double precision arithmetic. 

The numerical solution of the extended Orr-Sommeffeld equation (26) was obtained for 
a wide range of wavenumbers a,fov R equal to 5000, 11000 and 17000. At larger values 
of R the sizes of the arrays become prohibitively large, if good resolution is desired close 
to the wall. However, from (28), it is clear that if one is interested in wall modes it is not 
really necessary to compute at very high R. Indeed, it is seen from the results presented 
below that trends become nearly universal in the reported range of R. 

For all the three Reynolds numbers noted above, unstable eigenvalues were found for a 
wide range of a:. A typical set of results is shown in figures 7 and 8, for the point defined 
by (2.7, 5000) in the (a, R) plane. Here, Cr = 0.349 and a = 0.0044. The value of Cr 
indicates that this is a wall mode. (Note: u = 0.35 for y = 0.03, when R = 50(X).) Figure 7 
shows the eigenfunction (j>, normalised as <^ = 1 + 0/ at y = 1.5. Figure 8 shows plots of 
the r.m.s. values of u and v as estimated from the numerical solution. As seen from figure 8, 
r.m.s. values of u peaks at y+ 11, which agrees extremely well with experimental data 
(see Tennekes & Lumley 1972) and the results of direct numerical simulations (Spalart 
1988). Similar results were obtained for the production term, t/?(m)', where r/? = uv, 
which peaks at y''' 12.5, as shown in figure 9. (The inordinately high values of ordinates 

in figures 8 and 9 are no cause for concern, and follow from the normalisation (f) = I+ 0i 
at y = 1.5.) 

Figure 10 shows the growth rate of the organised disturbance, plotted against 

the wavenumber o;'^. One sees that the band of unstable wavenumbers is fairly large and 
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Figure 10. Graph of the growth rate, a+ct versus in inner variables scaling. 


virtually the same, in inner variables, for all the three values of R. Secondly, in keeping 
with (26), the growth rate curves seem to reach, asymptotically, a limiting curve. Most of 
the calculations reported are with X based on (23). However, one growth rate curve is also 
shown in the figure, with A. based on (22). As mentioned earlier, calculations using the 
latter expression of A, are stiffer and therefore not persisted with for higher R. Nevertheless, 
the growth rate curves show that the trend of the results, and the range of unstable a, are 
very comparable by these two widely different shapes of the A-curves. Thus two important 
conclusions drawn are as follows. First, the instability is governed by a large value of A 
(~ 9) in the inner region. Secondly, the exact shape of the A-curve is not important. 


5. Experiments 

The experiments were conducted in the 0.61 m x 0.61 m x 4.75 m test section, suction 
tunnel of the Gas Dynamics Laboratory of this department. Figure 11 shows a schematic 
drawing of the tunnel and the coordinate system used. The boundary-layer growing along 
the bottom wall of the test section was chosen for measurement. A boundary-layer trip 
(1.5 mm square rod) was placed at the start of the test section. Other details on the tunnel 
and the instruments used for measurements may be found in Sen et al (1993)and Sen & 
Veeravalli (1996a). 

Measurements were made at a nominal tunnel soeed of aDoroximatelv 8 m/s and at an x 
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the wall. The mean velocity profiles used in the numerical work and the best fit log-layer 
straight line from Purtell et a/ (1981) are also shown in figure 6. As mentioned before, 
since R is not very high, the inner log-law region does not extend over a very wide range 
in y'^. 

The friction velocity in the present experiments was estimated by several different 
methods. The best-fit log-layer straight line to various boundary layer experiments at 
comparable Reynolds numbers, reported in literature is (Purtell et al 1981); 

M+=5.0-1-5.62 log 10 (38) 

where, w"*" = and = y^Vi^/v = yR^fB. 

The value of for the present experiments may thus be estimated from both the slope 
and the intercept of the above log-law. From the slope one obtains = 0.352 m/s while 
from the intercept the value for is 0.353 m/s. The friction velocity may also be estimated 
from the Reynolds number Rq , using the skin friction data reported in literature (Purtell 
et al 1981). Using this method one obtains a value of 0.344 m/s. Finally, the Preston-tube 
measurement using the calibration of Patel (1965) yielded v-^ = 0.352 m/s. We see that 

Table 1. Boundary-layer parameters. 

8* = [ (1 — «/t/oo)d)'= 4.61 mm 
Jo 

pS - 

6=1 — u/[/oo)dy — 3.2Smm 

Jq 6^00 
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the different methods for estimating the friction velocity agree to within 2.5%. For the 
plots based on experiments in figures 5 and 6, has been taken as 0.352 m/s. 

The boundary layer thickness, defined for the experimental results as the location at 
which the mean velocity reaches 0.99f/oo, is 3.35 cm for this x location. Various parameters 
obtained from the mean velocity profile like the displacement thickness S* momentum 
thickness, 9 etc. are reported in table 1. The shape factor, H, falls well within the range 
reported in Purtell etal (19SI) and agrees to better than 1% with the data of Coles (1962), 
reported in Purtell et al (1981). 

Figure 12 shows compensated, one-dimensional, power spectra of the u' velocity fluc¬ 
tuations, /i:+£'J^j(/r+), plotted against /f+. Here k = 2nf/u, where / is the frequency. 



Figure 12. Graph of compensated power spectra k'^ E\\ (k'^) versus k'^, of the u' fluctu- 
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obtained growth rate curve of cx'^c^ versus k'^. Too much should not be read into the 
location of the ‘peak’ of this curve, as this curve only gives the growth rates according to 
linear stability theory. Nevertheless, the versus k'^ curve is contained fully within 
the collapsed experimental spectral curve, and within the energy-containing eddies region. 
This is yet another feature that keeps alive the possibility of some connection between the 
instability waves and actual turbulence. 


6. Conclusions 

An extended Orr-Sommerfeld equation has been derived to accurately describe the evo¬ 
lution of organised disturbances in the presence of background turbulence. This has been 
done based on a few novel features in the formulation. First it has been emphasized that the 
turbulent stresses in the problem can be suitably modelled based on the generalised eddy 
viscosity hypothesis (GEVH). This is a very significant observation to make, and this is 
the basis for legitimacy of the formulation of Reynolds & Hussain (1972), as well as of the 
present formulation. Second, after stipulating the validity of GEVH in the present problem, 
a significant and logical improvement has been made in modelling the Reynolds stress 
tensor, over Reynolds & Hussain (1972), based on the generalised anisotropic formulation 
of Pope (1975). This improvement is also unique in form. Thirdly, in the formulation, the 
anisotropy is properly accounted for in terms of a universal non-dimensional function A., 
which has been called the ‘anisotropy function’ herein. 

The numerical solution of the extended Orr-Sommerfeld equation (26) yielded unstable 
eigenvalues over a wide range of a. It has been found that instability is obtained by a large 
value of A(~ 9 or more) in the wall region. However, the calculations are not sensitive to 
the exact shape of the A,-curve. The mode of instability obtained is the wall-mode, and the 
eigenfunction and the r.m.s. distributions of u and v are found to be qualitatively similar 
to those obtained for Blasius flow. Moreover, the locations of the respective peaks of the u 
r.m.s. profile and the production term tr {u)', agree well with experimental values reported 
in literature, for actual turbulence. Thus the organised disturbances seem to mimic some 
of the features of actual turbulence. 

The band of unstable eigenvalues, when scaled with respect to inner variables, indicates 
that results approach universality with high R. Particularly, the growth rate curves for 
different R, in inner variables, viz. versus «+, approach a limiting curve for high 
R, as may be seen in figure 10. This is also a feature that correlates well with the inner 
scaling in actual turbulence. 

The experimental rig produces an adequately well-conditioned flat plate turbulent bound¬ 
ary layer flow to permit comparison with and validation of the numerical results. Detailed 
measurements of the one-dimensional power spectra for u' were taken at different in 
the buffer and log regions. These spectra were found to collapse well into a single curve. 
The limiting growth rate curve of the instabilities, plotted as a'^cf versus k^, is seen 
to be contained entirely within the power spectral curve of the u' fluctuations, and the 



entire range of unstable q:+ is contained within the energy containing eddies range. This is 
yet another feature that seems to support the possibility of some connection between the 
instabilities and actual turbulence. 

As a final word in conclusion, it may be stated that the present work keeps alive the 
question of a possible connection between stability theory and actual turbulence, in wall- 
bounded turbulent flows. 
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Abstract. The matrix condensation technique in conjunction with the sub¬ 
structuring principle has been previously used to model complex commercial 
automotive mufflers. Though complete 3-D finite element modelling is partially 
eliminated, it still requires the nodal-coordinate data and connectivity data of 
the elements forming the segment of each substructure. It also demands the 
connectivity of the degrees of freedom left after matrix condensation. Thus, the 
data preparation phase is tedious and cannot be claimed to be simple. More¬ 
over, the existing FEM code requires the user to have a knowledge of acoustic 
theory and finite element modelling. This paper describes a way to completely 
eliminate the data preparation phase. Given the geometric dimensions of the 
muffler as input, the muffler performance is obtained as the output. The results 
obtained are compared with analytical and experimental results for simple as 
well as extended-tube expansion chambers, with or without an offset. 

Keywords. Mufflers; duct acoustics; automatic mesh generation; finite ele¬ 
ment methods; noise control. 

1. Introduction 

Exhaust noise is a predominant source of noise in an automobile. Various types of muffler 
components are used to bring down the exhaust noise levels. Commercial exhaust mufflers 
have large cross-sections and complicated shapes. Therefore, 1-D analysis is not adequate 
for proper understanding of muffler performance. Thus, there is a need for a tool to perform 
a complete 3-D analysis of mufflers. The automotive industry is vibrant and fast-changing. 
Therefore, analysis tools used must be highly user-interactive and should be engineered 
to offer readymade solutions. The need for expertise on the part of the user should be 
eliminated so that industries need not rely upon experts to do muffler analysis. Thus, there 
is a need for a tool which will give muffler performance, given the geometric dimensions 
of the muffler. 
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While generally a one-dimensional (or, plane-wave) analysis (Munjal 1987a) is ade¬ 
quate, higher order modes do play an important role at higher frequencies in and around 
expansion chambers. Several methods have been proposed by the previous researchers 
for 3-D analysis of mufflers (El-Sharkawy & Nayfeh 1978; Ih & Lee 1985; Munjal 
1987b; Abom 1990; Sahasrabuddhe et al 1991). Analytical solutions proposed by re¬ 
searchers (El-Sharkawy & Nayfeh 1978; Ih & Lee 1985; Abom 1990) could be one 
possible solution. However, these methods cannot yield a completely automatic code 
as they require eigenfunctions which are not readily known for complicated shapes. 
Therefore, one has to resort to methods not involving eigenfunction expansions. 
Various numerical techniques proposed in the literature were the collocation method, 
finite difference method and the popular finite element method. The collocation method 
(Munjal 1987) is relatively simple in-so-far as the algebra and computational effort are 
concerned, but is deficient in dealing with fractional area ratios and complicated shapes. 
Therefore, the collocation method fails to be a solution for the problem on hand. 
Finite difference method is not suitable for very complicated and irregular geometries 
(Baumeister 1980). However, the finite element method is completely general as far as 
the complex geometry of the muffler and different properties of the medium are con¬ 
cerned, and hence can handle rigid/non-rigid walled chambers, and simple or extended tube 
chambers with inline/offset inlet/outlet pipes (El-Sharkawy & Nayfeh 1978; Zienkiewicz 
1979; Baumeister 1980; fli & Lee 1985; Munjal 1987; Abom 1990; Sahasrabuddhe et al 
1991). 

Finite element method was extended to acousto-structural analysis by Gladwell (1965). 
Later, Young & Crocker (1975) and Craggs (1976) used FEM for prediction of muffler 
performance. Young & Crocker (1976) used rectangular Hermitian elements with C\ con¬ 
tinuity while Craggs used the axisymmetric ring element obtained from the hexahedral 
element. Their FEM predictions compared well with experiments. Later, Christiansen &, 
Krenk (1988) used a recursive finite element technique for modelling pipes and ducts. 
They used an axisymmetric quadrilateral element and hence their analysis is restricted 
to 2-D configurations. Recently, Craggs (1989) analysed branched duct networks using 
FEM. A truly 3-D analysis of expansion chamber mufflers incorporating the recursive 
FEM technique of Christiansen & Krenk (1988) was carried out by Sahasrabudhe et al, 
(1991) who used the sub-structuring principle in conjunction with the matrix condensa¬ 
tion technique for modelling simple/extended tube expansion chamber mufflers. Though 
Sahasrabudhe et al{l99\) could analyse simple/ex tended tube expansion chambers, the 
FEM code was not user-interactive. The code requires the nodal coordinate data of the 
segments of the substructures, elemental connectivity data, and degrees of freedom left 
after condensation in addition to the dimensions of the muffler. Thus the data preparation 
phase is tedious and cumbersome. Moreover, the user should have expertise in acoustic 
theory, finite element modelhng etc. Thus the need for a truly 3-D analysis tool is not 
fulfilled. 

In this paper, the sub-structuring and the matrix condensation technique have been used 
to model muffler geometry. Given the dimensions of the muffler as the input, the 3-D 
finite element model is generated and muffler performance is obtained as the output. The 
results obtained using the automatic code are compared with the results available in the 



Figure 1. Acoustic system. 



2. Finite element formulation 

The acoustic field in a muffler configuration filled with stationary ideal fluid (figure 1) is 


given by 

[v2 + *2]p = 0. (1) 

with the following four types of boundary conditions; 

(1) ap/an = -jpoCD{p/Z), (2) 

when the normal acoustic impedance Z is known; 

(2) dp/dn - -jpocoun, (3) 

when the normal particle velocity is known as in the case of the constant velocity 
piston source; 

(3) dp/dn = 0, (4) 

at the rigid walls of the muffler; and 

(4) P = Pn, (5) 

when the constant pressure source is assumed. 

3 

Assuming a trial solution for the field variable p, at any point within an element as 

P = WeiPn}e, ( 6 ) 

where Ne is the basis function or the shape function or the interpolation function. Substi¬ 
tuting (6) into (1), we get the residual as 

Re = [V^Nl + k^Nj]Pne- 0) 

This error is orthogonalized with respect to the weighting functions to get 

( NeRe^v=^Q, ( 8 ) 

he 

which can be written as 

I Ne^^Nlpne^V+ ( k^NeN^Pne<^V 0. 


Ve ^ he 


(9) 



Then, following Munjal (1987a), we get 


[M],- 

k\S]e+^[D]e 

z J 

Prte — jPo^^Fei 

(10) 

[M]e = 

/ V iVg V nJ d V , (mass matrix), 

JVe 

(11) 

J 

I NgNj dV, (stiffness matrix), 

'Ve 

(12) 

[D]e = 

/ NgNg dS, (damping matrix), 

JSe 

(13) 

{F}e = 

/ UnNgdS, (forcing vector). 

JSe 

(14) 


Here, [S] denotes the stiffness matrix and the surfaces Se^, Se-, and denote the portion of 
the boundary in which normal acoustic impedance, the input source velocity and the rigid 
wall boundaries are known, respectively. Surface Se includes Se^, and 5^^. Equation 
(10) can be written after global assembly as, 

[K]{p} = {F} (15) 


where [F] is called the admittance matrix given by 

[*:i = = E 


(F) = = -jPoO> (17) 

i=l i=\ 

where / denotes the number of elements in the system. The contribution to forcing vector 
{F} comes only from the surface Si at the inlet of the muffler assuming constant velocity 
source and contribution to the damping matrix {D} comes from the only surface Sz at the 
outlet of the muffler (see figure 1). 


[M]e - kHsh + —[DL 


(16) 


3. Muffler performance using FEM 

Higher order modes or 3-D effects are generally limited to expansion chambers and a 
couple of diameters length of the inlet tube and outlet tube respectively. Therefore, a muffler 
composed of a volume chamber with inlet and outlet tubes (figure 2) can be simulated by a 
linear acoustic four-pole system (Munjal 1987a). Making use of this analogy, the acoustic 
state variables viz., acoustic pressure and volume velocity on the inlet side can be related 
to those on the outlet side through the four-pole parameters also called the transfer matrix. 
The transfer matrix relation can be written as 


'Pn 


[All 
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Figure 2. Engine exhaust system. 


where p and v are, respectively, acoustic pressure and volume velocity, subscripts n and 
1 stand for exhaust pipe and tail pipe (figure 2) respectively, and four-pole parameters are 
given by 


Aii = 

Ai2 = 
^21 = 
A 22 = 


Pi JD[ =0 

El' 

L 1^1 - pi=o 

Vn 

Pi Jj;i=0 
LtJ). p[=o 


(18) 

(19) 

( 20 ) 
( 21 ) 


The assembled set of dynamic equations ,(10), obtained through finite element formulation, 
are solved twice using two sets of boundary conditions to calculate four-pole parameters. 


(1) By imposing u = 0 (i.e.dp/dx = 0) at the outlet boundary, nodal points An and A 21 
are computed using (18) and (20). 

(2) By imposing p = 0 at the outlet boundary, nodal points A 12 and A 22 are computed 
using (19) and (21). 


Then the transmission loss (TL) value can be evaluated as, 

1 l.i , O ^12 , ^llPoCo , ^ 

TL = 20*log - * |Aii-h -1- - -I-A 22 I 

-2 PoCo Sn 


(22) 








where Sn and Si are the areas of the exhaust pipe and tail pipe respectively, and v is the 
acoustic volume velocity, 
po - density of the medium, 
and Co - sound speed. 

4. Selection of the finite element for muffler analysis 

In general, there is a wide-choice of finite elements available in the literature. However, 
the 20-noded isoparametric Brick element (Zienkiewicz 1979) has been used extensively 
throughout this work. The reasons for selecting Isoparametric 20-noded brick elements 
are as follows. 

(1) Since automotive mufflers involve curved geometry, 20-noded isoparametric elements 
with quadratic shape functions are found suitable. A large number of linear elements 
may be required to discretize the same geometry. Therefore, use of a higher-order 
element reduces memory requirement and solution time. 

(2) 20-noded isoparametric elements have been successfully used by Sahasrabudhe et al 
(1991) for finite element modelling of mufflers (Munj al 1987a) and the results obtained 
are satisfactory. 

5. Maximum length of a segment and element 

Investigation of the criteria for the maximum length of a typical segment is presented 
here. Ross (1981), based on his numerical experiments, reported that, in order to ensure 
convergence of finite element results up to a frequency of /-max, the maximum typical 
dimension of the finite element should be less than 0.2 times the minimum wavelength. 
He used 8-noded isoparametric elements in his mesh configuration. Later, Sahasrabudhe 
et aZ (1991) reported that, with the use of a quadratic 3-D element (20-noded isoparametric 
element), the maximum typical dimension of the finite element should be less than 0.52 
times the minimum wavelength. To verify this criterion, an analytical solution was derived 
for a simple square pipe and it was compared with the results using FEM. 

Details of the comparison are as follows: 

A square pipe with its boundary conditions is shown in figure 3. The dimensions of the 
pipe are such that only plane waves can propagate. Therefore, the governing equation of 
the problem is 

^ {D^^k^)p{z) = 0, (23) 


j~-*Z r--->A 


P(0)=1 


P(1)=0 


—>A 


SECTION A-A 


Figure 3. Square pipe model. 
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and the boundary conditions are; 

» @ z = 0, p{0) = 1, 

• @ z = l,p(l) = 0. 

The solution for (23) is 

p(z) = Cl ■ txp(-jkz) + C2 • &xp{3-jkz), (24) 

where Cl and C2 are constants and k is the wave number. The constants Cl and C2 can 
be found by applying the two boundary conditions. Thus, the solution is given by 

p{z) = s,mk{l + z) sin(^/), (25) 

where 

/ = (n .Co)/(4-/), for « = 1, 3, 5, 7 • • •. (26) 

Equation (25) is taken as a reference for validation of the FEM results. 

To verify the criteria for maximum element length, the following test runs were 
conducted. 

Data for the test run: 

(1) Frequency 1000 Hz, 

(2) length I = (n • Co)/(4 • /), n = 1,3, 5,..., 

(3) side of the square h < Co/(2 • /) 

The following procedure was adopted for conducting the test, and the steps involved are 
as follows; 





Figure 5. Comparison of FEM predictions (*) for a square pipe model of length 0.765 m, 
side 0.05 m and a segment length 0.225 times the wavelength with analytical solution (o) 
given by (25). 


• Select a suitable value of the pipe length using (26), 

• vary the segment length each time (greater than 0.52 times the minimum wavelength to 
start with), 

• the value of h (side of the square cross-section) is taken equal to the segment length so 
that each element has an aspect ratio of unity, 

• compare the results with the analytical solution. 

Analysis was conducted with a segment length equal to 0.7,0.54,0.375,0.225 times the 
minimum wavelength. The results for two cases are shown in figures 4 and 5. The solution 
obtained using 0.225 times the minimum wavelength (figure 5) tallied quite well with the 
analytical solution. From the above results it may be observed that the element length must 
be certainly less than 0.375 times the minimum wavelength. 

A limitation of the foregoing analysis is that it is conducted in the plane wave zone 
where the effect of higher order modes is not present. Hence, to check the validity of the 
above tests further, a complete 3-D analysis of a simple expansion chamber muffler was 
performed and the computed values of TL were compared with those obtained from an 
analytical solution. The tests are as below. 

A simple expansion chamber with a diameter ratio of 3 and length to diameter ratio of 
1.67 was analysed and the transmission loss thus obtained was compared with the results 
available in the literature (El-Sharkawy & Nayfeh 1978). The segment length was taken 
as 0.125 m. From figure 6, it may be observed that the results agree well up to a frequency 
of 865 Hz. This observation leads to the conclusion that the segment length should be 
less than or equal to 0.3 times the minimum wavelength of interest. Figure 7 shows the 
comparison between TL of a simple expansion chamber with the same dimension as in the 
previous case but differing in that the segment length taken is 0.062 m. The results agree 
quite well up to a non-dimensional frequency parameter kl 2 of 8 which corresponds to 





Figure 6. Comparison of TL of a coaxial expansion chamber with /i = I 3 = 0.375 m, 

di ~ d 2 = 0.05 m, I 2 = 0.2505 m, d 2 = 0.15rh, and with segment length 0.125 m, (-) 

FEM, (—) analytical (El-Sharkawy & Nayfeh 1978). 


a frequency of 1730 Hz. This observation confirms the previous conclusion. From these 
numerical experiments the following conclusions were derived. 

• A 20-noded isoparametric element can be used to model pipes with curved geometry. 

• Maximum typical dimension has to be less than 0.3 times the minimum wavelength 
when using the 20-noded isoparametric element and less than 0.2 times the minimum 
wavelength when using linear elements as given by Ross (1981). 



Figure 7. Comparison of TL of a coaxial expansion chamber with k = = 0.187 m, 

di = d^ = 0.05 m, h = 0.2505 m, di = 0.15 m, and with segment length 0.062 m, (-) 

FEM, (—) analytical (El-Sharkawy & Nayfeh 1978). 
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6. Modelling of simple expansion chambers 

A typical simple expansion chamber muffler is shown in figure 8. 

Basically it consists of three substructures (Munjal 1987a): 

• Inlet pipe (substructure 1); 

• expansion chamber (substructure 2); and 

• the outlet pipe (substructure 3). 

The task of a mesh generator is to discretize the volume of the muffler. The FEM code 
takes the length and diameter of the inlet pipe, the expansion chamber and the outlet 
pipe respectively as the inputs. The basis for the finite element modelling is the matrix 
condensation technique (Munjal 1987a). The construction of a 3-D finite element model 
proceeds as follows. 

Nodes are generated in a 2-D region as shown in figure 9 which uses only the inlet pipe 
diameter and the chamber diameter. The mesh shown is a cross-section of a typical segment 
of substructure 2. This 2-D mesh can be extruded to form a typical segment of substructure 
2 (Sahasrabudhe etal \99l). The 8-nodes on the periphery of the inner circle are extruded 
to form pipe segments for the inlet/outlet pipes. Number of segments of substructure 2 is 
given by m where it becomes m -f-1 for inlet/outlet pipes. The segment at the inlet section 
of substructure 1 and outlet section of substructure 3 are not condensed, as the gradient of 
pressure is computed using these segments. This is useful in the computation of four-pole 
parameters as discussed. 



Figure 9. Simple expansion chamber discretization. 
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Figure 10. Comparison of FEM results for a simple expansion chamber with d\ = 0.05 m, 
dijdi — 3, hjdi = 1.67 with the analytical solution of El-Sharkawy & Nayfeh (1978). (*) 
analytical, (o) FEM. 


Here, m is any positive integer such that the ratio 1/2"* is less than 0.3 times the 
inimum wavelength. One has to note that only nine elements are generated irrespective 
' the muffler dimension. This mesh is found to be sufficient to analyse mufflers of practical 
terest. Since the same number of elements are formed, the connectivity data for a typical 
jgment remain the same and are read from a data file. In the same manner, the information 
igarding the degrees of freedom left after matrix condensation is read from the data file, 
hus, the geometric dimensions alone are sufficient to obtain a 3-D finite element model 
f a simple expansion chamber muffler. Limitations of the code are as follows. 



Figure 11. Comparison of the TL for a simple expansion chamber with d\ = 0.0205 m. 
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Figure 12. Comparison of the FEM predictions for coaxial expansion chamber with d\ = 
0.0205 m, d 2 /d\ = 6, l 2 ld 2 = 1.35 with the experimental results of Ih & Lee (1985) (—) 
FEM, (o) experiments. 


• The circumference of the inlet/outlet pipe should be less than 1.2 times the minimum 
wavelength of interest. 

• The difference between the radius of inlet pipe and expansion chamber should not exceed 
0.3 times the minimum wavelength. 

These limitations may easily be overcome by using an adaptive mesh generator which 
will create more elements for a bigger dimension so that the maximum typical dimension 
of any of the elements does not exceed the prescribed limit. Thus the process of finite 
element modelling is fully automated without user intervention. 

To validate the modelling procedure discussed in the previous section, the results ob¬ 
tained using the completely automatic FEM code were compared with the experimental 
and analytical results available in the literature. They are as shown in the figures 10-12. The 
finite element results showed good agreement with the analytical solution (El-Sharkawy 
& Nayfeh 1978) throughout the frequency range of interest as can be seen in figure 10. 
From figures 11 and 12, it is seen that the FEM predictions are in good agreement with 
the experimental results as well. 


7. Modelling of extended tube expansion chambers 

A typical extended tube expansion chamber with five substructures is shown in figure 13. 
One may note that modelling of extended tube expansion chambers is very similar to that 
of simple expansion chambers. The difference lies in the additional task of modelling the 
two additional substructures (annular cavity introduced by the extended pipes). Figure 13 
shows the discretized cross-section of the annular cavity between the inlet pipe and the 
chamber. The nodes are formed in much the same way as the simple expansion chamber but 
the difference lies in the connectivity. In the armular cavity model, the element representing 
the cross-section of the inlet pipe is absent. It is to be noted here that there is no interaction 
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Figure 13. Discretized cross-sections of the annular cavity and the expansion chamber. 

between the substructures 1 and 2, or 4 and 5 as the walls of the pipe are assumed to be 
rigid. The interaction is only through substructure 3, The input to the FEM code will be the 
length and diameter of the five substructures. The FEM program can handle the following 
cases: 

• extended inlet and extended outlet, 

• simple inlet and extended outlet, 

• extended inlet and simple outlet. 

To validate the modelling procedure illustrated above, an extended tube expansion cham¬ 
ber was analysed and the results were compared with the experimental results of Ih & Lee 
(1985). The comparison is shown in figure 14. 
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Figure 15. Discretized cross-section of an offset inlet/outlet model. 

8. Modelling of offset inlet/outlet chambers 

In the case of offset pipe chambers, the automatic FEM code requires as the input the 
centre and radii of the pipes forming the substructures. The mesh generation routines are 
capable of handling any arbitrary position of the inlet and outlet pipes provided the two 
pipes do not intersect. Figure 15 shows the mesh generated for an angular position of 180° 
between the inlet and outlet pipes. 

Two-dimensional nodal data are converted into a segment model as shown in figure 
15, which after condensation is representative of the whole volume of the muffler. Thus 



Figure 16. Comparison of the FEM predictions (—) for an offset expzmsion chamber with 
the experimental results of Ih & Lee (1985). a = 6, ^ = 1.35; n = rs = 10.25 mm; 
Si = 82 = 031r2-, e = 180°. 














Figure 17. Comparison of FEM predictions (—) for an offset expansion chamber with 
the experimental results of Ih Lee (1985). a = 6 , ^ = 0.33; r\ = r^, = 10.25 mm; 
Si = <52 = 0.37r2; 6 = 180°. 


matrix condensation eliminates the task of constructing a complete 3-D model, which is 
quite cumbersome and consumes a lot of memory and CPU time. The nodal connectivity 
data for the segments of different substmctures are read from a data file. The number 
of segments in each of the substructures is calculated in such a way that the ratio //2"^ 
(segment length) is less than 0.3 times the minimum wavelength. Here m denotes the 
number of successive matrix condensations (Christiansen & Krenk 1988). Thus the process 
of modelling becomes fully automatic. 

Figures 16 and 17 corroborate the finite element results obtained in the form of TL of 
acoustically long as well as short offset inlet/outlet expansion chamber mufflers with the 
experimental results of Ih & Lee (1985). 



Figure 18. Comparison of the TL of the offset simple (—) and offset extended tube (—) 
expansion chambers with Id and leo = 0.33^2 for the extended tube, a = 6, = 2; Si = 0, 

82 = 0.63r2, 0=0°. 
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9. Modelling of offset extended inlet/outlet pipes 

The modelling of mufflers with inlet/outlet pipes offset and extended into the expansion 
chamber is done using the concepts discussed in the previous sections. Two more sub¬ 
structures are introduced by extending the inlet/outlet pipes. These cavities are modelled 
in the same way as the offset expansion chamber, the difference being the absence of the 
element corresponding to the pipe segment. While the element eleven corresponding to 
the inlet pipe is absent in the segment model of substructure two, element one is missing 
in the segment model of the substructure four (see figure 15). A typical muffler with axial 
inlet and radially offset outlet pipe was analysed with and without the extension of the 
inlet/outlet pipes and the results compared. The comparison is shown in figure 18. The 
resonances in the end chambers are found to improve the performance of the extended 
tube chambers as compared to their simple chamber counterpart. 

10. Integration of 3-D analysis with 1-D analysis 

Often, higher order modes that are generated in expansion chambers decay quickly in the 
adjoining inlet/outlet pipes, where only plane waves can propagate. Thus, the completely 
automatic 3-D analysis subroutines developed for an expansion chamber with small (about 
two diameters long) portions of the inlet pipe and outlet pipe can be used to generate a 
2x2 transfer matrix as described above. This can then be integrated with Munjals Transfer 
Matrix-based Muffler Program (TMMP) which, in fact, was the primary purpose of the 
investigation. Whenever expansion chambers appear as a part of the acoustic system, the 
3-D analysis routines are called upon to get the four-pole parameters connecting the inlet 
and outlet sections of the chamber. Thereafter, respective routines of the program are used 
to perform plane-wave analysis. The overall performance of the acoustic system can be 
obtained by simply multiplying the successive transfer matrices (Munjal 1987a). As an 
illustration, consider an acoustic system which consists of a simple expansion chamber and 
a concentric tube resonator in series as shown in figure 19. The overall performance of the 
system is calculated and plotted in figure 20. It may be noted that 1-D analysis results show 
good agreement up to the cut-on of the first azimuthal mode, but show discrepancies in the 
high frequency range. This observation clearly illustrates the inadequacy of plane-wave 
analysis above the cut-on frequencies of higher order modes. But the CPU time for each 
frequency step in the case of a 3-D FEM program is at least 100 times that of a plane-wave 
analysis program. Therefore, one has to resort to 3-D analysis only when the higher order 
mode effects predominate within the frequency range of interest. Thus it is possible to 
analyse any combination of simple/extended tube chambers in series using the TMMP 
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Figure 18. Comparison of FEM predictions (—) with TMMP (—) for the acoustic system 
shown in figure 19. 


with 3-D analysis routines being restricted to chambers whenever the user feels that 3-D 
analysis is inevitable. 


11. Concluding remarks 

• Maximum typical dimension of a finite element should be less than 0.3 times the min¬ 
imum wavelength of interest while using quadratic shape functions and 0.2 times the 
minimum wavelength for linear shape functions, as suggested by Ross (1981). 

• A completely automatic FEM program which takes dimensions of the muffler as the input 
and four pole parameters as the output has been developed for analysing simple/extended 
tube expansion chambers with coaxial or offset location of the inlet/outlet pipes and 
integrated with Munjal’s TMMP software. This work thus fulfills the need for a user- 
friendly tool for 3-D analysis. It is to be noted that the need for an expert with knowledge 
of acoustic theory and finite element modelling is totally eliminated. 
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Abstract. This paper presents an overview of various exact analytic and ap¬ 
proximate numerical methods for the solution of radiative heat transfer prob¬ 
lems in participating media. Review of each method is followed by its strengths 
and limitations. Importance of radiative heat transfer analysis and difficulties 
in the solution of radiative transfer problems have been emphasized. 
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1. Introduction 

The rapidly depleting fossil fuel reserves, the growing ecological imbalance due to com¬ 
bustion products, the requirements of the space age etc., have necessitated the designing 
of energy efficient systems, and this in turn requires accurate prediction of heat transfer 
rates in these systems. Of the three modes of heat transfer, conductive and convective heat 
transfer rates are generally linearly proportional to temperature differences, whereas ra¬ 
diative heat transfer rates are mostly proportional to differences in the fourth powers of the 
absolute temperatures. Thus, in high-temperature applications, determination of radiative 
heat transfer rate becomes an important and necessary consideration. Consideration of 
thermal radiation becomes paramount in the design of high temperature energy conversion 
devices (boilers, engines, rocket nozzles etc.), in industrial processes (heating furnaces, in¬ 
cinerators), in nuclear reactors (such as in the Sun, fusion reactors, nuclear bombs), during 
atmospheric reentry of space vehicles, in materials processing (glass, metals), in utiliza¬ 
tion of materials (metallic screens, ceramic foams), in analyses of industrial and building 
fires etc. Thus, this wide range of potential applications provides strong motivation for 
continued research in the field of radiative heat transfer. 
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Especially in combustion applications, either for power production or propulsion, con¬ 
sideration of thermal radiation becomes extremely important because combustion phe¬ 
nomenon as well as heat transfer rates are greatly influenced by radiation. For example, in 
large-scale coal combustion chambers (Manickavasagam & Menguc 1993) and in fluidized 
bed combustion (Tien 1988), radiative heat transfer accounts for as much as 90% and 40% 
respectively of the total heat transfer rates. Further, accurate determination of radiative 
heat transfer is essential in the design of burners, furnaces, and internal combustion (IC) 
engines. In these systems, minimization of emission of the oxides of nitrogen, NOjt, is 
one of the primary concerns, because NO^ formation is highly sensitive to the changes in 
temperature (Kuo 1986). 

Although in high temperature applications radiation remains the dominant mode of heat 
transfer, for evaluation of total heat transfer rate, conduction and convection also have 
to be considered. In the solution of any high temperature flow problem, which involves 
conduction, convection, and radiation (conjugate problem), apart from the usual difficul¬ 
ties encountered due to the large number of partial differential equations required for the 
solution of conductive and convective parts of the problem, a greater complexity in prob¬ 
lem formulation as well as much higher computational cost become imminent due to the 
presence of an extra radiative heat flux term in the general equation for energy conservation. 

In full simulation of the heat transfer process in a typical high-temperature energy 
conversion device, say an internal combustion engine, one has to consider three momentum 
equations, at least two turbulent equations, equations for species, and one general equation 
for energy conservation incorporating all the three modes of heat transfer. The general 
equation for energy conservation for moving compressible fluid is given as 


/ar 

pe„(-+v.VT 


- V • ( -/rVr -h ) -pW . V -I- -b Q'". 

conduction radiation 

»- _ ■ _ - ^ ^ -- - _ ^ 

divergence of total heat flux 


( 1 ) 


In the general energy conservation equation (1), the term on the left hand side represents 
the contribution to energy transfer owing to convection, whereas the terms on the right hand 
side respectively represent contributions of conduction, radiation, pressure work, viscous 
dissipation, and internal generation. 

Even though in (1), radiative heat flux is part of the total heat flux term, it itself is the 
most difficult part to evaluate, particularly in the presence of participating media. The 
radiative part has to be determined from the equation of conservation of radiative energy, 
which is an integro-differential equation, and is given as 

V • q/? = - G(rx)], (2) 

where the blackbody intensity, in,, in (2), is the Planck’s intensity function which depends 
explicitly upon the local temperature, T, of the medium and kx is the local linear spectral 
absorption coefficient of the medium. The incident radiation, G, in the above equation, is 
given by 



where in the above equation, the intensity of radiation at any optical depth, tx, in an 
absorbing, emitting, and scattering medium is given by 

i(rx, 9, 0) = i(0, 9, 0) exp(-rx) + f S(tI, 9, (f)) exp(-(rA - t{)) dz^. 

J 0 

(4) 

Equation (4) is the integral form of the radiative transfer equation (RTE). It describes 
the intensity of radiation at some optical depth xx in the medium in terms of the intensities 
reaching that position from all other locations x-^ in the medium and from the medium 
boundary at = 0. It is obtained by integrating the differential form of the RTE (given 
below) from the point of origin of the ray, = 0 at the wall, to a point x^ = xx inside the 
enclosure. 

di{xx., 9, 0) 

— + Krx, (t>) = S(Tx, 0). (5) 

dr 

The source function in (5) is given by 

S{xx,9,(t)) = {l-(D)ixb{tx)+{-^\ ( f i{xx,9' 

\4:;r / J(p=o Je^o 

xp((9',c()') (^,0))sin0'd0 d0', (6) 

where the single scattering albedo, cox = crxs/iicx + and axs is the scattering coeffi¬ 
cient, and p((9 (9, 0)) is called the anisotropic scattering phase function which 

represents the probability distribution of the energy scattered into the direction (9, 0) from 
the rays traversing the elemental volume from all other directions, (9 ,0 ). 

It should be noted that the RTE given in (4) is a third-order integral equation in intensity, 
i{xx,9,<p). For its evaluation, the integral over the source function must be carried out over 
the optical coordinate, xx- But the source function itself is an unknown function of intensity 
(see (6)) which has to be integrated over two directions, 9 and 0. Furthermore, usually the 
temperature T which is related to the blackbody intensity term, ixb, in the source function, 
(6), is not known and it must be found in conjunction with the overall conservation of 
energy, (1), thus further complicating the situation. 

Therefore, due to the presence of radiation, the general equation for energy conservation, 
(1), turns out to be an integro-differential equation which poses enormous difficulties in its 
solution. Further, complexities due to radiation crop up due to the fact that radiative heat 
transfer deals with seven independent variables, viz., three space coordinates, the frequency 
of radiation, two coordinates describing the direction of photon travel, and time, whereas 
conduction and convection involve a maximum of four independent variables, viz., three 
space coordinates and time. Thus, the combination of the above complexities makes the 
complete solution of the problem, in the real sense, a formidable task even with the recent 
advent of supercomputers. Thus, there has always been the quest for the development of 
efficient and cost-effective methods, for predicting radiative heat transfer rates in practical 
energy conversion systems. 

A number of numerical techniques exist in the literature for the determination of radia¬ 
tive heat transfer rates in various engineering geometries with participating media. Many 




new methods. As real problems are almost always of a conjugate nature, according to 
Howell (1983) and Yang (1986), any practical numerical radiation method should have the 
following features. 

(1) Capability to handle multi-dimensional and complex enclosure geometry. 

(2) Good accuracy under all conditions: non-scattering or anisotropic scattering, gray or 
non-gray, isothermal or non-isothermal participating media. 

(3) Ease of application. 

(4) Ease of generalization: flexibility in choosing the different order of approximations. 

(5) Availability of the accurate intensity field and the integrated quantity. 

(6) Computational compatibility with conductive and/or convective heat transfer codes. 

(7) Low computation cost. 

Therefore, the objective of the research related to the development of radiative heat transfer 
models for determination of heat transfer rates should be aimed at achieving the above 
criteria. 

The purpose of this review is to acquaint the reader with different approaches in prac¬ 
tice presently for solving problems of radiative transport, applicable to heat transfer, in 
absorbing, emitting, and scattering media. There exists a large body of literature relevant 
to radiative transport in participating media, and it is not feasible to mention them all in any 
paper. Also, in these days of numerous journals and allied publications, it is quite possible 
that some of the major works might have been inadvertently overlooked. Further, any review 
process is rather an arbitrary process because the authors have to decide what to include, 
what to omit, and where to start and end. These points, in one sense, reflects the authors’ 
biases, and frankly, the present article is not an exception. In the following sections, refer¬ 
ence has only been made to the significant publications related to various methods. A brief 
introduction of each method has been provided along with its strengths and limitations. 

2. Exact analytical method 

Analytical methods can be applied to radiative heat transfer problems only in the case 
of highly idealized situations, like problems with simple geometry and homogeneous 
participating media having spectrally independent radiative properties. In such situations, 
exact analytical solutions can be obtained for the integro-differential radiative transfer 
equation, which is used for determining the variation of intensity in the participating 
media. RTE applied to one-dimensional plane parallel media have abundant applications 
in various fields, such as atmospheric sciences, neutron transport etc., where their exact 
analytical solutions have received much attention (Viskanta & Menguc 1987). 

However, exact analytical solutions of RTE for radiative heat transfer in 1-D enclosures 
and for simple situations are now the subject matter of text books on thermal radiation 
(see, e.g., Ozisik 1973; Siegel & Howell 1981; Modest 1993). For one-dimensional cylin¬ 
drical geometry, Heaslet & Warming (1966) have given an exact solution for the case of an 


isotropically scattering medium enclosed in an isothermal black cylinder. Azad & Modest 
(1981) have provided exact analytical solutions for gray, linearly anisotropically scatter¬ 
ing media confined in one-dimensional cylindrical enclosures with arbitrary but specified 
temperature distribution. Loyalka (1969) has also solved a similar kind of problem and his 
method is also considered exact. 

Cheng (1968) has used a rigorous approach to solve the RTE in a two-dimensional 
rectangular enclosure with absorbing-emitting-isotropic scattering media. Crosbie and 
his coworkers (Crosbie & Koewing 1979; Crosbie & Dougherty 1980, 1985; Crosbie & 
Schrenker 1982, 1983, 1985; Crosbie & Farrell 1984; Crosbie & Davidson 1985), have 
provided exact formulations of RTE for the case of absorbing-emitting-anisotropically 
scattering media for two- and three-dimensional rectangular and cylindrical enclosures. 

For complex radiation problems (for example, problems with inhomogeneous participat¬ 
ing media, spectrally dependent radiative properties etc.), the exact analytical solutions of 
the integro-differential RTE are exceedingly difficult, and exphcit solutions are impossible 
for all but very simple situations. Therefore, for such complicated problems, approximate 
numerical methods are necessary. These methods have been reviewed in the following 
sections. Although, exact analytical methods can be adopted only for some simple ge¬ 
ometries with many simplifying assumptions, they can be useful in giving quantitative 
indications for more difficult situations, and serve as benchmarks against which many 
other approximate numerical methods can be checked. 

3. Zonal method 

The zonal method, usually known as Hottel’s zone method, is one of the most widely used 
methods for calculating radiation heat transfer in practical engineering systems. It was first 
developed by Hottel & Cohen (1958) for an absorbing-emitting and non-scattering gray 
gas with constant absorption coefficients. Hottel & Sarofim (1967), Noble (1975), Smith 
et al (1985) and Nelson et al (1986) extended it to deal with non-constant and non-gray 
absorption coefficients as well as isotropically scattering media. 

In this method, the bounding walls and the enclosure volume are divided into a number 
of zonal elements. Each element is assumed to have a uniform temperature distribution and 
radiative properties. Then the direct exchange areas, between each surface and volume 

elements are evaluated and the total exchange areas, which include the view factors for the 
direct exchange as well as those due to multiple reflections, SiSj, are detennined. Using 
these exchange areas, the zonal method expresses the radiant balance of each element in 
the form: 

N 

Qi = QE.i-Y^Sl^Ej. (7) 

j 

Here the net radiant heat flow, Qi, from element i, is the difference between emitted 
radiation, Qej, and the total radiation received directly and indirectly through the multiple 
reflection and scattering from all the N elements in the domain. The term Ej is the total 
blackbody emitted flux, crTj^, and 5/5), the total exchange area from element j to i, is 
a function of the direct exchange area, J[sJ, and the radiative properties of the medium. 


Once the temperature field of the medium is specified, using the matrix inversion technique, 
algebraic equation (7) can be solved. 

As originally formulated tlie zonal method has some inherent limitations, such as the 
treatment of non-gray, temperature-dependent radiative properties (Hottel & Cohen 1958). 
If the radiative properties are temperature-dependent, the exchange areas have to be re¬ 
calculated over and over again, as the temperature field is determined iteratively. For com¬ 
plex geometries, application of the zonal method becomes difficult, as numerous exchange 
factors between zones must be evaluated and stored in the computer memory. Further, it 
is often difficult to couple the equations of the zonal method with momentum and gen¬ 
eral equations for energy conservation, (1), which are usually solved using finite element, 
finite difference or finite volume techniques. This is mainly because of the different sizes 
of the required control volumes used in solving radiative equations and momentum and 
general equations for energy conservation. A coarser grid is desirable for the zonal method. 
This method can be completely prohibitive if the finer grid scheme as used by the finite 
difference equations is adopted (Viskanta & Menguc 1987). 

Seeing the capability of the zonal method to yield nearly analytical solution, researchers 
have been trying to enhance its applicability by coupling it with momentum and energy 
equations. Steward & Tennankore (1979) have coupled the zonal method with finite differ¬ 
ence equations in modelling a combustor by adapting two different grid schemes; one for 
the radiation pait and the other for the flow and the temperature field. Smith et al (1986) 
and Sistino (1982) have coupled the zonal method with momentum and energy equations 
to solve the combined radiative and convective heat transfer in an absorbing-emitting 
and isotropically scattering medium flowing through a cylindrical duct. More recently, 
Yuen & Takara (1994) have modified the zonal method to deal with the anisotropic scat¬ 
tering case, but its application to conjugate problems with anisotropic scattering is still 
awaited. 


4. Monte Carlo method 

The Monte Carlo method (MCM) is a purely statistical sampling technique. It can exactly 
simulate almost all important physical processes without much difficulty. In this method, 
the numerical treatment of the mathematical formulation for incorporation of non-gray, 
temperature-dependent gas properties and anisotropic scattering situations is easy. Also, 
the usual difficulties encountered in complex geometries can be circumvented easily. It is 
because of these advantages that the MCM has been applied to solve many radiative heat 
transfer problems. 

In this method, radiation is assumed to be composed of discrete energy bundles or 
photons. The computation consists of following the probable path of a finite number of 
photons, which obey the physical restraints imposed by the RTF, until their final absorption 
in the system. The path of tlie photon is determined at each point of emission, absorption- 
emission, scattering or reflection by a random choice from a set of possible paths. A 
sufficiently large sample of photons are followed to get statistically stable results. 

This method has been used extensively in atmospheric and neutron transport studies. 
In heat transfer studies, it was first applied to one-dimensional radiative heat transfer 
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problems in early 1960s by Fleck (1961), Howell & Perlmutter (1964a, 1964b) and 
Perlmutter & Howell (1964). Radiative heat transfer problems of increasing complexities, 
such as non-gray, anisotropic scattering media confined in various types of enclosures, that 
have been investigated by this method have appeared in the works of Taniguchi (1967, 
1969), Stockham & Love (1968), Avery etal (1969), Gupta et al (1983), Abed & Sacadura 
(1983), Yang et al (1983), Tiwari & Liu (1992), Farmer & Howell (1994), and Liu & 
Tiwari (1994). 

Though this method is highly accurate, it is computationally very expensive. To make 
the Monte Carlo method more attractive for use in conjugate problems, Kobiyama et al 
(1979) and Kobiyama (1986) have studied ways of reducing the computational expense. 
In certain test problems, they could reduce the computational time by up to 15%. For 
coupled conduction-radiation heat transfer in semi-transparent media. Abed & Sacadura 
(1983) have used the finite-difference Monte Carlo method. Liu & Tiwari (1994) have 
solved conjugate problems for the case of non-gray participating media using narrow band 
model. Farmer «& Howell (1994) have obtained Monte Carlo solutions for radiative heat 
transfer in a three-dimensional enclosure with anisotropically scattering inhomogeneous 
media at constant temperature. Most recently, Mishra & Blank (1995) have provided Monte 
Carlo solutions as a benchmark for two-dimensional enclosures with absorbing-emitting- 
isotropically scattering media at radiative equilibrium, 

Monte Carlo solutions are always very close to the exact analytical solutions, but due 
to the statistical nature of the method, they fall within a narrow band around the exact 
analytical solution. Thus due to its accuracy, the method has extensive application to large 
varieties of engineering problems, and hence it has been concluded that the MCM should 
be used as a benchmark for comparison of other methods (Viskanta & Menguc 1987; 
Mishra & Blank 1995). However, like other numerical methods, MCM also has some 
disadvantages. Its large appetite for computer time, inability to match the required grid 
size needed for concomitant computation of conduction and/or convection parts of the 
conjugate problems, and statistical fluctuation of results, limits its use to benchmarking of 
only purely radiation problems. 


5. Multi-flux method 

In multi-flux method, to simplify the governing RTF, angular (directional) dependence 
of intensity is decoupled with its spatial dependence. Intensity is assumed to be uniform 
over a given interval of the solid angle (Lockwood & Shah 1976; Khalil et al 1982). The 
integro-differential form of RTE is then averaged over each of the finite angular regions 
to result in a set of flux equations. 

Depending on the number of discrete directions over which the solid angle is discretized, 
different multi-flux methods result. If the solid angle is divided uniformly into two, four 
or six divisions, the result is two-flux, four-flux and six-flux methods, respectively. The 



flux approximations in a three-dimensional cylindrical enclosure containing an absorbing, 
emitting, and scattering medium. 

In general, the accuracy of the flux approximation depends on the choice of the solid 
angle subdivisions. If there is no intersection between the adjacent subdivisions, more 
accurate results are expected (Abramzon & Lisin 1984). For rectangular enclosures, the 
same has also been observed by Siddal & Seluck (1979). If the distribution of radiation 
intensity is assumed for each subdivision, the general equation given by Abramzon & 
Lisin (1984) can be simplified and solved simultaneously. If the fluxes in each subdivision 
are assumed constant, a simpler six-flux model can be obtained from the general flux 
equations. For an absorbing, emitting and scattering medium in a cylindrical enclosure, 
Spalding has suggested a simpler six-flux model. A different variation of the six-flux model 
has also been proposed by Chu & Churchill (1960). Although, Chu and Chirchill’s six-flux 
model is for a one-dimensional plane parallel medium, it is possible to modify it for multi¬ 
dimensional enclosures. Starting from this six-flux model, for axi-symmetric cylindrical 
enclosures, Varma (1979) has obtained a four-flux model. 

There are mainly three objections to the multi-flux approximations, for practical en¬ 
gineering problems. Firstly, there could be no coupling between the axial and the radial 
fluxes, which makes the equations physically unrealistic. Secondly, the approximation of 
the intensity distribution from which the flux equations are obtained is arbitrary. Lastly, 
the model equations cannot approximate highly forward scattering correctly, even though, 
theoretically it is feasible (Viskanta & Menguc 1987). 

6. Moment method 

In moment method, radiation intensity is expressed by a truncated series representation in 
terms of the products of angular and spatial functions as 

N 

y> z? 0) = Ao -f -|- rf Afi^y (8) 

n=\ 

where A 5 are functions of locations only and r), and jj, are direction cosines in x, y, z 
directions, respectively (Viskanta & Menguc 1987). Although, this equation is written in 
Cartesian coordinates, it can be given for any orthogonal system. 

This method is the generalization of the Milne-Eddington approximation to the higher 
order. The Milne-Eddington approximation involves an approximation to the angular de¬ 
pendence of intensity such that the integro-differential RTE is transformed into an ordinary 
partial differential equation (Modest 1993). 

7. Spherical harmonics method 

The spherical harmonics method or P// approximation uses a series approximation to 
the angularly dependent intensity field to replace the integro-differential RTE with a set 



of partial differential equations (Menguc & Viskanta 1985, 1986; Ou & Kuo-Nan 1982; 
Ratzel & Howell 1983). The intensity field is expandedin terms of the spherical harmonics, 
where N terms are retained for the Pj\/ approximation. 

Like most other approximations, this method was also originally developed for the study 
of radiative heat transfer in the atmosphere. Latter on, it was modified for solution of neutron 
transport problems (Davison 1958), and extensively used for one-dimensional heat transfer 
problems (Kourgnaoff 1952; Chandrasekhar 1960; Ozisik 1973; Cess & Sparrow 1978; 
Siegel & Howell 1981). Bayazitoglu & Higenyi (1979) have given formulations and solved 
radiative heat transfer problems for non-scattering Cartesian, cylindrical and spherical 
media using the first-order P] and third-order P 3 spherical harmonics approximations. 
For an isotropically scattering cylindrical medium, Higeny & Bayazitoglu (1980) and 
Ahluwalia & Im (1981) have used the P\ approximation. Ratzel & Howell (1983) have used 
Pi and P 3 approximations to solve isotropically scattering two-dimensional rectangular 
medium problem. Menguc & Viskanta (1985, 1986) have given the general formulation 
of the Pi and P 3 approximations for an absorbing, emitting, and anisotropically scattering 
medium in two-dimensional, finite cylindrical as well as three-dimensional rectangular 
enclosures. Liu et al (1992) have analysed the effect of anisotropic scattering on radiative 
heat transfer in a two-dimensional rectangular enclosure using Pi approximation and 
have compared their results with higher order discrete ordinate method employed by Kim 
& Lee (1988). It has been shown by Krook (1955) that the moment method and the 
spherical harmonics method are closely related and completely equivalent (using spherical 
harmonics, which are functions of direction cosines, and exploiting their orthogonality 
properties). 

The spherical harmonics method or Pm approximation is one of the most tedious and 
cumbersome of the radiative transfer approximations (Viskanta & Menguc 1987). The 
formulation of the spherical harmonics method can be very complicated even for the low 
order P 3 method, in multi-dimensional geometries. Errors can easily result during the 
long derivations (Menguc & Viskanta 1985, 1986; Ou & Kuo-Nan 1982). Generalization 
to higher order approximations is not realistic with this method. The Pi approximation is 
very accurate, if the optical dimension of the medium is large (i.e., greater than 2) however, 
it yields inaccurate results for thinner media (Viskanta & Menguc 1987), especially near 
boundaries. The Pi approximation becomes less reliable, if the radiation field is anisotropic, 
i.e., there are large temperature and/or particle concentration gradients in the medium. The 
P 3 approximation, however, can yield accurate results for an optical dimension as small 
as 0.5 and for anisotropic radiation fields at the expense of additional computational effort 
(Menguc & Viskanta 1986; Ratzel & Howell 1983). 


8. Discrete ordinate method 

The discrete ordinate method or the Sm method was first proposed by Chandrasekhar 
(1960) in his work on stellar and atmospheric radiation. Originally, this method received 
little attention in the heat transfer community. Again, like other methods, this method 
was first systematically applied to problems in neutron transport theory, notably by Lee 
(1962) and Lathrop (1966). In the beginning, there had been some unoptimized attempts to 



apply this method to one-dimensional, planar thermal radiation problems by Love & Grosh 
(1965), Hsia & Love (1967), Hottel et al (1968) and Roux & Smith (1974). However, this 
method has been applied only very recently to, and optimized for, general radiative heat 
transfer problems through the pioneering work of Truelove (1987, 1988), Fiveland (1984, 
1987,1988), Kim & Lee (1988), Cheong & Song (1995), and Vallion et al (1995, 1996). 

This method is obtained by discretizing the entire solid angle (^2 = 47r) using a finite 
number of ordinate directions and corresponding weight factors. The N discrete values 
of direction cosines Cn, l^n always satisfy the identity -f- -f = 1. The RTE 
is written for each ordinate direction and the integral terms are replaced by a quadrature 
summed over each ordinate direction. The ordinate directions are usually the Gaussian or 
Lobatto quadrature points, although different ordinates and weights have been proposed 
(Fiveland 1987; O’Dell & Alcouffe 1987). The integral over the angles is then replaced by a 
quadrature sum. The resulting set of angularly discretized RTE can be spatially discretized 
by using finite difference method. 

This approximation is also not flawless. One of the major drawbacks of this method is the 
presence of ray-effects, which lead to anomalies in scalar flux distribution (Lathrop 1968; 
Lewis & Miller 1984). These effects are caused by the inability of the discretized intensity 
distribution to represent the actual continuous distribution fully, especially when radiant 
energy conservation is not insisted upon. Consequently, radiation can be lost if it does not 
fall into one of the discrete ordinate directions. Unless an infinite number of directions are 
used, ray effects cannot be totally eliminated. The ray effects are especially pronounced 
if there are localized radiation sources in the medium and scattering is less important in 
comparison to absorption. As the single scattering albedo increases, the radiation field be¬ 
comes more isotropic and the ray effects become less noticeable. However, with increasing 
single scattering albedo and/or optical thickness of the medium, the convergence rate may 
become very slow (Viskanta & Menguc 1987). Also, the computer time and memoiy re¬ 
quirements of this method are relatively high as compared to Pn approximation and zonal 
method (Kim 1990). 

Another drawback of the discrete ordinate method is the difficulty in extending its 
application beyond Cartesian geometry (Pomraning 1973). This is evident even in the 
simple cylindrical system where formulation is complicated by additional restrictions by 
the presence of curved surfaces (Carlson & Lathrop 1968; Hyde & Tmelove 1977). 

It should be noted that the basic principle of this method and the multi-flux method are 
the same, but the difference is in the evaluation of the integral. The flux method uses a 
trapezoidal rule for the angular integral, whereas the discrete ordinate method uses a higher 
order Gaussian type quadrature formula, which is more accurate (Khalil et al 1982). In 
the limit of a large number of discrete directions, both methods will converge to the same 
result. 


9. Discrete transfer method 

This method, proposed by Shah (1979) and Lockwood & Shah (1981), combines the virtues 
of the zonal, Monte Carlo, and discrete-ordinate methods. In this method, calculation of ra¬ 
diative heat transfer involves the tracing of representative rays from one surface to another 



through the domain of interest. Every bounding control surface in the computational do¬ 
main is the source of a number of rays chosen to divide the hemisphere about a grid point 
in a number of equal solid angles. The intensity distribution along each ray is calculated 
by solving a discretization of the RTE. The RTE along a ray, neglecting scattering can be 
expressed in the form (Shah 1979) 


di kgT^ 

~ -Ki A - 
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(9) 


where i is the intensity of radiation, s is the direction along the ray, jc is the absorption 
coefficient, and a is the Stephan-Boltzmann constant. For any representative ray, the 
intensity distribution can be calculated by assuming each control volume as homogeneous. 
Under this assumption, above equation is integrated to give the recurrence relation. 
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where in and in+i are the intensity of the ray on entry and exit, respectively, and As is the 
distance travelled in the control volume. Tc is the temperature in the centre of the control 
volume. Therefore, given the initial intensity at a point on the emitting surface, the change 
in intensity along the ray can be calculated using (10). The initial intensity is specified by 
taking the walls to be Lambert surfaces. 

For a gray wall, the emitted intensity is dependent on the incident flux q ~ 
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( 11 ) 


where iy^(6,(j)) is the incident intensity. In discrete transfer method, the incident flux 
integral is replaced by a numerical quadrature 

' q~ = '^i~{6,(p)cos9 sinOAOA^. (12) 


The values of incident intensity, i~, at a point on a wall are calculated by tracing rays from 
the point and backtracking from the walls intercepted, applying (10) through each control 
volume. Hence, for a grey-walled enclosure, the coupling between the incident flux and 
emitted intensity make the discrete transfer method a guess and correct procedure in which 
an estimate of the incident flux distribution is iteratively improved. 

Lockwood Sc Shah (1981) showed that very accurate results could be obtained with this 
method in one- and two-dimensional geometries by increasing the number of rays. More 
recently, da Graca & Fontes (1993) have solved radiative heat transfer problems in two- and 
three-dimensional rectangular enclosures with absorbing, emitting, and isotropic scattering 
media. Their comparison of the results, with other well established methods, such as the 
zone method, the discrete ordinate method and the method, shows reasonably good 
agreements. 

More recently, Murthy & Choudhary (1992) have employed discrete transfer method 
to convection-radiation problems in general complex enclosures having non-scattering 
participating media. Most recently. Cumber (1995) has suggested some improvements to 
the discrete transfer method and his suggestions have shown some improvements in the 



result. However, the improvement in Cumber’s (1995) results over that of Shah’s (1979) 
does not seem significant. 

Although, this method is claimed to be capable of accounting for anisotropic scattering 
in the medium, so far no results have been reported or compared against other bench¬ 
mark methods for scattering media in multi-dimensional enclosures. The results for a 
one-dimensional scattering medium did not show the same level of agreement with the 
benchmark results as did the non-scattering medium predictions. 


10. Other methods 

There are a few methods whose developments have not yet progressed to a level whereby 
definitive evaluations can be made. However, for the sake of completeness of the review, 
they are included here. 

10.1 Heat ray method 

This approach to the solution of radiative heat transfer equation for multi-dimensional 
enclosures has been presented by Taniguchi et al (1986) for absorbing-emitting media. 
This method is based on the Bouguer’s law and yields radiant energy absorption distribution 
in non-isothermal enclosures containing combustion gases. Comparisons of the predictions 
based on this method with other results show that the method is more accurate and less 
time-consuming than both zonal and Monte Carlo methods, if the radiative properties such 
as the absorption coefficient and wall emissivities are constant (Taniguchi et al 1986). 

10.2 Discretized intensity method 

This method has been proposed by Shih & Chen (1983) in solving a two-dimensional 
conduction-radiation problem in absorbing-emitting participating media. It was later used 
for solving a conjugate convection-radiation recirculating flows problem in a 2-D enclosure 
(Shih & Ren 1985). Unlike conventional approaches, in which radiative fluxes are assumed 
to be one-dimensional, this method bypasses the unnecessary use of integro-differential 
energy transport equation and accounts for two-dimensionality of the radiative fluxes. In 
its formulation, the conditions at an individual control volume have an explicit effect on 
the four immediately adjacent control volumes. This method has been shown to be only 
valid for relatively optically thick participating media (Blank 1992). 

10.3 Finite element approach 

The finite element method has been applied by Razzaque et al (1983, 1984), Sokman 
& Razzaque (1987), Chung & Kim (1984) and Fernandes et al (1981), to the combined¬ 
mode conduction and/or convection problems with radiation, including scattering and both 
known temperature and known flux boundary conditions. The finite element method, in 
principle, gives exact solutions within the errors introduced by the numerical scheme itself, 
i.e., no approximation is required in the formulation of a problem. According to Howell 



Tan (1989) has proposed and applied the product integration method to radiative heat trans¬ 
fer problems. This method, while analogous to zone method or the finite element method, 
requires evaluation of only three- or two-dimensional integrals for three-dimensional sys¬ 
tems. This reduction in the dimension of the required integration, significantly reduces the 
required computation time. This gain is achieved because for a given degree of approxi¬ 
mation of the temperatures (i.e., constant, linear, binomial etc.) within each of the n finite 
volume or surface elements, the number of calculations in the zonal and finite element 
method increases as while in the product integration method, it increases as n (Howell 
1988). Although Tan (1989) has incorporated linear anisotropic scattering effect in the 
analysis of a two-dimensional absorbing, emitting, and scattering media, the formulation 
is extremely complex. 

10.5 Finite volume approach 

The finite volume approach for radiative heat transfer problems proposed by Raithby & 
Schneider (1988) and Raithby & Chui (1990) shares the same philosophy as the finite 
volume technique currently used in predicting fluid flow and convective heat transfer 
problems. Therefore, it can potentially use the same non-orthogonal computational grid as 
the fluid flow solver. Conceptually, the method can handle non-uniform radiative properties, 
anisotropic scattering, and accounts for spectral effects. The application of this method 
to the radiative heat transfer problems in one- and two-dimensional Cartesian (Raithby 
& Chui 1990) and axi-symmetric (Chui 1990) enclosures, with absorbing, emitting, and 
anisotropic scattering media, has shown good comparisons with the existing methods. 

10.6 Collapsed dimension method 

The collapsed dimension method is the latest development which has been proposed by 
Blank (1994) and Blank & Mishra (1995,1996). Unlike the rest of the radiation models, it 
views radiative phenomena quite differently. This method collapses the 3-D radiative in¬ 
formation of the problem into its 2-D solution plane. To accomplish this task, the collapsed 
dimension method makes use of the effective intensity. Each effective intensity contains 
in itself the information of the conglomeration of all the conventional intensities, focused 
on the same point and contained in a plane which is both inclined at the same angle as 
the effective intensity and also normal to the solution plane. This novel procedure thus 
eliminates the use of the solid angles and reduces both the complexity and computational 
expense involved in the solution of such problems. This method has been shown to give 
near analytical results for all ranges of optical thicknesses, including optically very thin 
cases where most of the existing methods fail to perform well. Since this method views 
all radiative phenomena in the 2-D solution plane, the computational time taken by this 
method for solution of any particular problem will be the least. For real life conjugate 
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problems, this method seems to be one of the best. However, this method has not yet been 
tested for such complicated situations. 


11. Conclusions 

Based on the above discussions, it is evident that no single radiation model can be regarded 
as the best model for all types of problems. Each method has its own merits and demer¬ 
its. Thus, to ascertain whether a particular method is appropriate for a given problem, 
comparisons have to be made against the benchmark results obtained either from exper¬ 
iments or exact analytical solutions. Although, the Monte Carlo method’s applicability 
to conjugate problems is very restricted, for the majority problems dealing purely with 
radiative transport, its solutions serve as benchmark results. The zonal method solution is 
also considered exact and has been successfully applied to a large class of problems. The 
Pn method and the discrete ordinate method, even being approximate methods, have the 
most diverse applications, in all types of simple and complex problems (conjugate prob¬ 
lems). However, these methods do not perform well at low optical thickness situations. 
The discrete transfer method has been found to give good results for different types of 
enclosures with simple boundary conditions and simple participating media conditions 
but for complicated situations it has not been tested. Finite element and finite volume ap¬ 
proaches to radiative transport problems have the added advantage in conjugate problems 
in which most of the methods have difficulties in matching the grid sizes of radiative part 
of the problem with conductive and/or convective part. The collapsed dimension method 
has been found to yield very accurate results for the problems tested so far and is believed 
to take least time compared to the existing methods. Unlike, the P^ and discrete ordi¬ 
nate methods, this method works extremely well at very low optical thickness conditions. 
Moreover, this method has not been tested so far for a large class of problems and thus 
no conclusive remarks can be made about complex situations. In conclusion, depending 
upon the computing facility and desired accuracy of the result etc. for any problem, it is 
the user’s discretion to adopt a particular method. 


List of symbols 

Cu specific heat at constant volume, [W/m^ kg K]; 
E emissive power, [W/m^ ]; 

Eb blackbody emissive power, [W/m^]; 

G incident radiant energy, [W/m^]; 

i intensity, [W/m^]; 

ib Planck’s blackbody intensity, [W/m^]; 

k thermal conductivity, [W/mK]; 

N number of volume elements in zonal method; 

p phase function; 

0 energy flow rate, [Wl; 


s geometric path, [m]; 

s[sj direct exchange areas in zonal method, [m^]; 

Si Sj total exchange areas in zonal method, [m^]; 

T temperature, [K]; 

t time, [s]; 

V velocity vector, [m/s]; 

X, y, z Cartesian coordinates, [m]; 
rj direction cosine, [-]; 

0 zenith angle, [rad]; 

T optical thickness, optical coordinate, [-]; 

xi enclosure optical thickness, [-]; 

K absorption coefficient, [m“ ^ ]; 

jjL direction cosine, [-]; 

jji dynamic viscosity, [kg/m s]; 

s; direction cosine, [-]; 

o Stephan-Boltzmann constant = 5,670 x 10~^ W/m^K"^; 

Os scattering coefficient, [m“^]; 

(p azimuthal angle, [rad]; 

d> dissipation function, [J/kg m^]; 

(jo single scattering albedo, [-]; 

p density, [g/m^]. 

Subscripts 

a absorption; 

b blackbody; 

5 scattering; 

X wavelength dependent. 

Superscripts 

+ positive direction; 

— negative direction; 

' dummy variables. 
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